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SECTION  1 


OVERVIEW  AND  SUMMARY 


1.1  INTRODUCTION 

This  report  describes  research  performed  during  the  second  year  of  a 
three  year  investigation  of  wave  propagation  and  related  effects  in  multiphase 
porous  media.  The  research  effort  is  comprised  of  both  a  theoretical  and  an 
experimental  program  which  compliment  and  support  one  another.  This  combined 
approach  has  proved  valuable  in  identifying  and  analyzing  phenomena  governing 
the  propagation  of  high  intensity  stress  waves  in  saturated  porous  media  and 
other  related  phenomena  including  liquefaction  and  post-liquefaction  processes 
such  as  consolidation.  The  report  is  organized  into  four  chapters  (Sections  2 
through  5)  covering  progress  in  theoretical  analysis  and  implementation  of 
analytic  techniques  into  computational  algorithms,  and  three  chapters 
(Sections  6  through  8)  describing  experimental  results  and  analyses.  An  over¬ 
view  and  summary  of  these  seven  chapters  is  provided  in  the  following  subsec¬ 
tion. 


1.2  CHAPTER  OVERVIEW  AND  SUMMARY 

Sections  2  and  3  contain  a  theoretical  and  numerical  two-phase  analysis 
of  undrained  uniaxial  strain  and  hydrostatic  loadings  of  porous  saturated 
materials.  Uniaxial  strain  loadings  closely  duplicate  exp 1 os ion- induced 
airblast  loadings  on  the  ground  surface  as  well  as  direct  coupled  loadings 
along  a  spherically  divergent  wavefront  at  early  times.  Uniaxial  loadings  are 
thus  the  most  generally  applicable  loadings  for  the  study  of  explosive 
effects.  Section  2  discusses  the  advantage  of  the  multiphase  analysis  from 
both  the  theoretical  and  experimental  standpoints,  and  develops  the  theoreti¬ 
cal  constitutive  equations  for  undrained  loading  of  a  fully  coupled  elasto- 
plastic  skeleton  with  compressible  grains  and  pore  water.  Specific  examples 
of  hydrostatic  and  uniaxial  loadings  are  shown  which  demonstrate  the  interre¬ 
lationships  between  effective  stress,  pore  pressure  and  total  stress  and  which 
demonstrate  liquefaction  early  in  the  unload  cycle  resulting  from  the  hystere- 
tic  nature  of  the  material  skeleton. 
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Section  3  describes  the  development  of  the  numerical  code  NKOCP  which 
uses  the  theoretical  equations  of  Section  2  in  incremental  form  to  model  the 
uniaxial  and  hydrostatic  undrained  response  of  actual  soils  and  rock.  The 
code  incorporates  nonlinear  elastic  response  of  the  solid  grains,  nonlinear 
compressibility  of  the  pore  water  and  nonlinear  plastic  response  of  the  porous 
skeleton  into  a  fully  coupled  two  phase  material  model.  Skeleton  response 
characteristics  from  laboratory  tests  can  be  input  directly  into  the  code  to 
insure  accurate  depiction  of  the  actual  material  response.  The  use  of  NKOCP 
is  described  and  a  listing  of  the  code  and  an  example  problem  are  provided. 

Sections  4  and  5  contain  theoretical  and  numerical  treatments  of  wave  pro¬ 
pagation  and  damping  in  saturated  porous  media.  In  Section  4,  a  closed  form 
solution  for  wave  propagation  velocity  and  damping  in  fully  saturated  porous 
media  is  derived  for  a  fully  coupled  model  with  compressible  solid  grains  and 
pore  water.  This  solution  demonstrates  existence  of  two  types  of  compression 
waves,  termed  waves  of  the  first  and  second  kinds.  It  also  provides  a  means 
of  benchmarking  and  verifying  multiphase  code  calculations  and  a  means  of 
planning  and  guiding  a  further  investigation  of  these  two  types  of  waves  using 
our  finite  element  two  phase  codes. 

The  general  theoretical  solution  for  the  wavespeeds  and  damping  are 
incorporated  into  the  numerical  code  TWAVE,  described  in  Section  4  and  listed 
in  APPENDIX  C.  TWAVE  is  used  in  Section  5  in  a  parametric  study  of  the 
influence  of  excitation  frequency  and  variations  in  material  properties  on 
propegation  velocity  and  damping.  Compressional  wave  velocity  for  waves  of 
the  fVsv  k'Jm*  is  shown  to  vary  as  a  function  of  the  frequency-permeability 
product,  with  a  zone  where  wavespeed  transitions  from  a  lower  bound  value  to  a 
higher  bound  value  with  increasing  values  of  the  product.  Damping  is  seen  to 
be  a  maximum  where  the  rate  of  change  in  wavespeed  is  greatest.  Waves  of  the 
second  kind  also  show  a  transition  in  wavespeed  from  near  zero  at  low  values 
of  the  frequency-permeability  product  to  an  upper  bound  value  at  higher  values 
of  the  product.  The  results  of  this  parametric  study  will  be  used  in  the 
coming  year  to  identify  and  study  phenomena  controlling  damping  and  wavespeed 
for  both  types  of  waves. 
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Descriptions  of  the  experimental  apparatus  and  results  are  contained  in 
Sections  6  through  8.  Section  6  shows  the  results  of  undrained  uniaxial 
strain  loading  of  saturated  limestone  and  uncemented  silt-sand-gravel  mix¬ 
tures.  Both  materials  show  liquefaction  early  in  the  unload  portion  of  the 
cycle  predicted  by  the  theoretical  model  developed  in  Sections  2  and  3. 
Comparisons  are  shown  between  the  laboratory  data  and  calculations  of  the 
undrained  response  made  with  NKOCP  using  completely  independent  measures  of 
skeleton  response  and  grain  and  pore  water  compressibility.  There  is 
excellent  overall  correlation  with  the  laboratory  data  and  the  computed 
undrained  material  response. 

Section  7  describes  an  experimental  study  of  grain  damage  during  various 
types  of  loading  to  delineate  similarities  and  differences  in  response  of  gra¬ 
nular  materials  on  the  microscopic  level  between  various  shear  loadings  and 
purely  compressional  loadings.  Grain  damage  was  found  to  vary  substantially 
with  variation  in  mean  stress,  but  has  a  relatively  minor  dependence  on  the 
level  of  shear  stress.  A  method  for  quantitatively  describing  both  the  amount 
and  severity  of  grain  damage  is  developed.  The  stress  strain  response  of  all 
materials  tested  in  this  study  is  included  for  use  in  future  analysis  in  which 
the  plastic  work  will  be  correlated  with  the  grain  damage. 

Section  8  describes  experimental  apparatus  and  results  being  used  to 
investigate  and  characterize  fluid  friction  in  saturated  porous  media.  Our 
generalized  model  for  fluid  friction  is  described,  which  includes  frictional 
energy  dissipation  in  the  laminar,  transitional  and  turbulent  flow  regimes. 

In  this  phase  of  the  study,  Biot's  theoretical  work,  which  was  restricted  to 
the  laminar  flow  regime,  is  verified  experimentally  in  both  circular  and  flat 
ducts.  Empirical  coefficients  are  obtained  which  model  the  fluid  fric¬ 
tion  behavior  in  the  transitional  and  turbulent  flow  regimes  as  well. 
Preliminary  tests  are  also  documented  showing  the  influence  of  fluid  inertia 
and  Biot's  acceleration  sensitive  dissipation  term  on  the  pore  pressure  gra¬ 
dient  during  nonsteady  state  flow. 
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SECTION  2 


THEORETICAL  TREATMENT  OF  UNDRAINED  UNIAXIAL  AND  HYDROSTATIC  LOADINGS 

2.1  INTRODUCTION 

In  order  to  model  the  undrained  response  of  saturated  porous  soils  and 
rocks  to  high  intensity  dynamic  loadings,  it  is  often  advantageous  to  combine 
the  skeleton  properties  with  the  properties  of  the  pore  fluid  and  solid  grains 
using  two  phase  constitutive  relationships.  Specific  advantages  of  this 
approach  are: 

•  it  avoids  the  experimental  difficulties  of  achieving  100% 
saturation  and  avoids  the  inaccuracies  in  measuring  effec¬ 
tive  stress  by  differencing  total  stress  and  pore  pressure 
data.  The  latter  inaccuracies  arise  because  effective 
stresses  are  often  very  small  relative  to  the  total  stresses 
and  pore  pressures. 

•  it  permits  computation  of  reliable  material  response  charac¬ 
teristics  in  high  pressure  regimes  which  are  extremely  dif¬ 
ficult  to  achieve  in  the  laboratory  but  which  are  important 
in  the  study  of  explosive  effects  or  other  high  energy  appli¬ 
cations. 

Another  application  of  the  two  phase  theoretical  model  is  in  the 
interpretation  and  understanding  of  laboratory  data  from  undrained  load-unload 
testing  of  saturated  soils  and  rocks.  As  is  described  in  Section  6,  we 
liquefied  both  porous  soil  and  rock  in  undrained  (biaxial  strain  tests.  The 
two  phase  theoretical  model  was  very  useful  in  interpretation  of  these  test 
data,  and  allowed  us  to  identify  some  unexpected  aspects  of  pore  pressure  and 
effective  stress  response.  This  new  insight  into  two  phase  response  resulted 
in  revision  of  some  of  our  earlier  work  based  on  a  more  elementary  two  phase 
model.  It  also  increases  our  understanding  of  blast  liquefaction  phenomena 
which  have  been  observed  in  large  scale  field  tests  (Blouin  and  Shinn,  1983). 
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2.2  IDEALIZEO  RESPONSE  OF  FULLY-COUPLED  SATURATED  MATERIALS  TO  UNDRAINED 

HYDROSTATIC  AND  UNIAXIAL  STRAIN  LOADINGS. 

Fragaszy  and  Voss  (1986)  observed  liquefaction  of  saturated  sands  during 
the  unloading  portion  of  undrained  hydrostatic  tests.  Kim,  Blouin  and  Timisn 
(1986)  reported  prel irainary  test  results  showing  liquefaction  of  saturated 
soils  and  rocks  during  the  unloading  portion  of  undrained  uniaxial  strain 
tests.  Their  results  are  updated  and  completed  in  Section  6  of  this  report. 
The  liquefaction  phenomena  during  a  single  load-unload  cycle  in  which  stresses 
are  either  completely  or  predominantly  compressive  rather  than  shear  can  be 
explained  by  two  phase  analysis  which  considers  the  inelastic,  hysteretic 
nature  of  the  soil  or  rock  skeleton  (Rischbieter  et  al.,  1977,  Blouin  and  Kim 
1983).  While  the  overall  liquefaction  response  could  be  explained  by  simple 
decoupled  two  phase  models,  subtleties  in  the  pore  water  and  skeleton  response 
were  incorrectly  described  by  these  models.  In  order  to  more  fully  understand 
the  liquefaction  data  presented  in  Section  6,  we  extended  this  earlier  work 
through  the  use  of  a  fully  coupled  two  phase  model. 


2.3  HYDROSTATIC  RESPONSE 


The  hydrostatic  and  uniaxial  strain  response  of  saturated  porous 
materials  is  described  by  Blouin  and  Kim  (1984).  The  three  governing 
equations  describing  the  hydrostatic  response  are  the  effective  stress  law 
expressed  as 

P  *  P*  ♦  R  (2-1) 

the  relationship  between  volume  strain  and  pore  pressure  and  effective  stress 
given  by 


Kg 


(2-2) 


and  the  bulk  undrained  stress-strain  relationship  given  by 


where 


p  =  total  pressure 
p*  =  effective  pressure 
jt  =  pore  pressure 
sv  =  volume  strain 
Ks  =  bulk  modulus  of  skeleton 
Kg  =  bulk  modulus  of  solid  grains 

Kf  =  fully  coupled  undrained  bulk  modulus  of  saturated  two 
phase  material 


Kf  is  given  by 


Kf  *  ^m  +  Ks 


KmKs 


Km^s 


k  4.  v  KmKs 
*m  +  *s  _ 


Kg2  -  KmKs 


(2-4) 


where  Kro  is  the  bulk  modulus  of  the  solid  grain-pore  water  mixture  given  by 

KaKw 

K”  “  K«  ♦  n?Kg  -  Kh)  <2"6) 

with 

Kw  <*  bulk  modulus  of  pore  water 
n  ■  porosity 

Eliminating  the  effective  pressure  from  Equation  2-1  and  2-2  yields 

p«7r^l-^j  +  Ksev  (2-6) 

Substitution  of  Equation  2-6  into  2-3  yields  the  pore  pressure  as  a  function 
of  the  volumetric  strain  as 
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(2-7) 


Eliminating  the  pore  pressure  from  Equations  2-1  and  2-2  gives 


P  =  P’  ^  +  Kgev 


(2-8) 


Substitution  of  Equation  2-8  into  Equation  2-3  gives  the  effective  pressure  as 
a  function  of  the  volumetric  strain  as 


(2-9) 


For  purposes  of  an  idealized  simple  model,  the  pore  water  and  solid 
grains  may  be  assumed  linear-elastic,  and  the  skeleton  response  can  be 
described  by  a  bilinear  approximation  as  shown  in  Figure  2.1  with  a  loading 
modulus  KSf  and  an  unloading  modulus  Ksu.  For  a  given  maximum  volume  strain, 
evmax*  th®  corresponding  maximum  effective  pressure  is  calculated  from 
Equation  2-9  as 


P’max 


(Kg  -  Kf|) 


evmax 


(2-10) 


where  Kfi  is  computed  from  Equation  2-4  using  K$  «  Ksj.  During  unloading 
the  effective  pressure  is  expressed  as 


(evmax  "  €v) 


(2-11) 


P‘ 


P'max  “ 


(Kg  -  Kfy) 


If,  during  undrained  unloading,  the  effective  pressure  drops  to  zero,  a  state 
of  liquefaction  occurs.  For  a  load-unload  cycle  beginning  with  no  pressure  on 
the  saturated  material,  the  volume  strain  at  the  onset  of  liquefaction,  evliq* 
can  be  calculated  from  Equation  2-11  with  p'  set  to  zero  as 


£->)  . 

evliq  =  evmax  "  (Kg  -  Kfu)  P  max 


(2-12) 


Substituting  Equation  2-10  into  2-12  yields  the  strain  at  the  onset  of 
liquefaction  as  a  function  of  the  peak  strain 


evliq 


(K«  - Kf<1  (h  - ') ' 

(Kg  -  Kfu)(^  -  l) 


e 


vmax 


(2-13) 


To  illustrate  the  behavior  of  a  saturated  material  in  hydrostatic 
undrained  loading  a  simple  numerical  example  is  presented  in  Figure  2.2.  The 
following  values  which  roughly  approximate  the  properties  of  a  porous 
limestone,  were  used  in  calculating  the  total  pressure,  effective  pressure  and 
pore  pressure  load-unload  cycles  shown  in  this  Figure. 

Kg  •  5.0  k  106  psi 
KS|  «  0.5  x  108  psi 
Ksu  *  2.0  x  108  psi 
Ku  •  0.3  X  10s  psi 
n  »  0.1 

From  Equations  2-5  and  2-4  the  following  values  of  the  mixture  modulus,  K*, 
and  the  undreined  loading  and  unloading  modulus  of  the  bulk  mixture,  Kfj  and 
Kf u,  were  computed  • 


9 


K,,,  =  1.948  x  106  psi 
Kff  =  2.142  x  106  psi 
Kfu  =  2.831  x  106  psi 

Assuming  a  maximum  volume  strain  of  3%,  values  of  peak  total  pressure,  pore 
pressure  and  effective  pressure  from  Equations  2-3,  2-7  and  2-9  respectively 
were  computed. 


Pmax  *  64,260  psi 
*max  ■  54,720  psi 
p'max  =  9 >540  psi 

Ouring  unloading,  because  of  the  hysteretic  nature  of  the  skeleton,  the 
material  liquefies  at  a  strain  computed  from  Equation  2-13  of 

£vliq  =  2.341% 

Using  the  values  comput'd  a»ove,  the  load-unload  response  of  the 
idealized  saturated  po**ous  material  is  plotted  in  Figure  2.2  in  terms  of  total 
pressure,  pore  pressure  and  effective  pressure  versus  volume  strain.  The 
first  significant  feature  of  this  plot  is  *Hat  liquefaction  occurs  at  a  strain 
of  2.34%  as  the  effective  pres'ure  between  the  grains  falls  to  zero. 
Liquefaction  as  a  result  of  the  hysteretic  skeleton  is  clearly  illustrated 
here.  The  skeleton  hysteresis  also  result?  in  hysteresis  in  the  bulk  material 
and  apparent  hysteresis  in  the  pore  water.  Note  that  the  apparent  hysteresis 
in  the  pore  water  is  negative,  i.e.  there  is  an  apparent  net  energy  gain  in 
the  pore  water  during  the  load-unload  cycle.  This  apparent  energy  gain  is  due 
to  rapid  expansion  of  the  solid  grains  during  unloading  which  results  from  the 
rapid  drop  in  effective  pressure  in  the  skeleton.  As  the  solid  grains  expend 
they  cause  a  less  rapid  decrease  in  pore  pressure  down  to  the  point  of 
liquefaction.  Once  liquefaction  is.  reached,  the  pore  pressure  unloads  at  a 
rate  governed  by  the  bulk  modulus  of  the  solid  grsin-rore  water  mixture,  Kg. 

The  hysteresis  in  the  total  pressure  load-unload  cycle  represents  the 
actual  energy  lost  by  the  bulk  material  during  the  cycle.  The  bulk  material 
loads  up  along  a  slope  given  by  the  loading  bulk  modulus,  Kfj*,  and  unloads  down 


a  slope  given  by  the  unloading  bulk  modulus,  KfU,  until  the  point  of  liquefac¬ 
tion  is  reached.  From  there,  the  material  is  a  dense  fluid  and  unloads  down 
the  slope  given  by  the  mixture  modulus,  Km.  Note  that  the  area  (energy  dissi¬ 
pation)  within  this  cycle  is  considerably  less  than  the  area  within  the  effec¬ 
tive  pressure  cycle.  The  area  within  the  effective  pressure  cycle  represents 
both  the  nonrecoverable  energy  dissipated  by  the  skeleton  as  heat,  and  the 
recoverable  energy  associated  with  the  elastic  compression  and  expansion  of 
the  solid  grains.  This  recoverable  portion  of  the  energy  equals  the  area 
within  the  pore  pressure  cycle;  thus,  the  total  area  within  the  effective 
pressure  cycle  equals  the  sum  of  the  areas  within  the  total  pressure  and  pore 
pressure  cycles. 

2.4  UNIAXIAL  STRAIN  RESPONSE 

The  uniaxial  loading  is  representative  of  the  strain  path  under  airblast 
loadings  as  well  as  the  early  time  strain  path  from  spherical  explosive 
charges.  For  uniaxial  loadings  Blouin  and  Kim  (1984)  present  derivations 
which  parallel  those  of  the  previous  subsection  for  hydrostatic  loadings.  The 
three  governing  equations  for  the  uniaxial  response  include  the  effective 
stress  law. 


a a  ■  oa'  +  n  (2-14) 

the  relationship  between  volume  strain,  effective  stress  and  pore  pressure, 


* 

«v  “  ♦ 


Ks 

KgMs 


n 


and  the  constrained  undrained  stress-strain  relationship. 


(2-15) 


where 


oa  Hf€v 


(2-16) 


cra  *  total  axial  stress 
cra'  *  effective  axial  stress 
Hs  *  constrained  modulus  of  skeleton 

Nf  a  fully  coupled  undrained  constrained  modulus  of  saturated 
two  phase  material 


11 


Mf  is  given  by 


Mf  *  Kf  +  Ms  -  Ks 


(2-17) 


Eliminating  the  effective  axial  stress  from  Equations  2-14  and  2-15 
yields 


oa  =  it 


i1  ‘  Kg  )  +  Ms£v 


(2-18) 


Substitution  of  Equation  2-18  into  Equation  2-16  gives 


(2-19) 


Note  that  substitution  of  Equation  2-17  into  2-19  yields  Equation  2-7.  Thus, 
the  pore  pressure  response  in  uniaxial  loading  is  identical  to  that  in 
hydrostatic  loading. 

Eliminating  the  pore  pressure  from  Equations  2-14  and  2-15  gives 


ffa  *  V 


(2-20) 


Substitution  of  Equation  2-20  into  2-16  gives 


Because  radial  stresses  are  generally  less  than  axial  stresses  in  :» 
uniaxial  loading,  an  additional  expression  for  effective  radial  stress  as  a 
function  of  volume  strain  must  be  developed.  Xn  their  development  of 


12 


governing  equations  for  undrained  uniaxial  loading,  Blouin  and  Kim  (1984)  give 
effective  radial  stress,  ar',  as 


ffr' 


W 


(1  -  2V)  ^S 
(1  -  y)  Kg 


n 


(2-22) 


where  the  coefficient  of  horizontal  earth  pressure,  K0#  is  given  by 


(2-23) 


and  v  is  Poisson's  ratio. 

Substitution  of  Equations  2-19  and  2-21  into  2-22  yields 


(l-v)(Kg-Ks) 


[yMsKg  +  MgKs  (1  -  2v)  -  MfKs  (1  -  V)]  ev  (2-24) 


In  order  to  develop  a  simple  idealized  model  for  uniaxial  loading  we 
assume  the  same  bilinear  model  described  in  Figure  2-1  for  the  hydrostatic 
loading  with  the  constrained  skeleton  modulus  given  by  the  elastic  rela¬ 
tionship 


Ms  •  "I'B-V)  “s  (2-26) 

Assuming  Poisson's  ratio  to  be  a  constant  during  both  loading  and  unloading, 
the  appropriate  constrained  loading  and  unloading  moduli,  M$f  and  Nsu  respec¬ 
tively,  can  fee  computed  from  the  hydrostatic  moduli  using  Equation  2-25. 

For  a  given  maximum  axial  or  volumetric  strain,  eVmax*  corresponding 
maximum  effective  axial  stress,  0'araax,  the  maximum  effective  radial  stress, 
o'rmax’  the  maximum  pore  pressure,  nnax,  and  the  maximum  total  axial  stress, 
0anax  are  computed  using  Equations  2-21.  2-24,  2-19  and  2-14  respectively  and 
the  appropriate  values  of  loading  moduli  Mfj.>  Msj  and  Ksj.  Similarly,  the 
effective  axial  and  effective  radial  stresses  during  unloading,  the  pore 
pressure  dur4 unloading  and  the  total  axial  stress  during  unloading  can  be 
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expressed  using  the  above  respective  equations  and  the  corresponding  unloading 
moduli,  MfU,  Msu  and  Ksu. 

In  order  to  compute  the  point  of  liquefaction  during  the  uniaxial 
loading,  Cy/Hq/  the  skeletal  conditions  at  the  point  of  liquefaction  must  be 
defined.  During  undrained  unloading  the  effective  radial  stress  and  effective 
axial  stress  do  not  generally  tend  to  reach  zero  at  the  same  strain.  However, 
as  soon  as  any  one  of  the  principal  effective  stresses  drops  to  zero  the  ske¬ 
leton  enters  in  an  unconfined  state  and  is  assumed  to  be  unable  to  support 
skeletal  stresses  in  the  other  principal  stress  orientations.  Thus,  we  define 
the  point  of  liquefaction  in  undrained  uniaxial  unloading  as  the  point  where 
any  principal  effective  stress  first  drops  to  zero.  Generally,  the  radial 
effective  stress  will  reach  zero  before  the  axial  stress.  Using  a  parallel 
development  to  that  for  Equation  2-13  for  the  hydrostatic  loading,  the  volume 
strain  at  the  onset  of  liquefaction  during  uniaxial  unloading  can  be  computed 
by; 

(1)  expressing  the  peak  effective  radial  stress  in  terms  of 
peak  strain  and  the  appropriate  loading  moduli  according 
to  Equation  2-24 ; 

(2)  expressing  effective  radial  stress  during  unloading  in 
terms  of  the  volume  strain  Uvmax  ’  cv)  and 
appropriate  unloading  moduli  using  the  expression 

for  o'rmax  *ro,n  U)  above  in  the  same  manner  used  to 
obtain  Equation  2-11; 

(3)  Solving  the  above  expression  for  the  volume  strain  at 
liquefaction,  ev^q,  by  setting  the  effective  radial 
stress  to  zero,  thus  obtaining: 


£v1iq  • 


(Kg-Ksu)  (vHslKq+Hs<KS|{l-2v)-Hf  fKsl(l-t>) ) 

(kg-KgJ  J  {  MsU^g^SuKsU  (  1“2V  )  -NfyK$y  ( 1”V)T 


6vmax 


(2-26) 
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For  a  fluid,  where  radial  stress  equals  axial  stress  (no  shear  stresses), 
Poisson's  ratio  is  0.5  and  Ms  equals  Ks  according  to  Equation  2-25;  Equation 
2-26  then  degenerates  to  Equation  2-13  for  hydrostatic  load-ings. 

To  illustrate  the  behavior  of  a  saturated  material  in  uniaxial  undrained 
loading  a  numerical  example  is  shown  in  Figure  2.3  using  the  same  material 
parameters  as  used  in  the  hydrostatic  loading  of  Figure  2.2.  The  only  addi¬ 
tional  parameter  needed  is  Poisson's  ratio  which  is  taken  as 

v  =  0.2 

From  Equations  2-25  and  2-17  the  constrained  loading  end  unloading  moduli  of 
the  skeleton  and  bulk  saturated  material  are  computed  as 

Msj  «  1.0  x  106  psi 
Nsu  a  4.0  x  106  psi 
Nf|  «  2.642  x  106  psi 
Hfu  «  4.831  x  106  psi 

Assuming  a  peak  axial  volume  strain  of  3%,  values  of  peak  total  axial  stress, 
port  pressure,  and  effective  axial  and  radial  stress  were  computed  from 
Equations  2-16,  2-19,  2-21  and  2-24  respectively  as 

Oamax  ■  79,260  psi 
irttax  8  64,730  psi 
ff'amax  *  24,630  psi 
ff'rmax  *  2,030  psi 

The  hysteretic  nature  of  the  skeleton  results  in  liquefaction  once  the  radial 
effective  stress  drops  to  zero  at  a  strain  given  by  Equation  2-26  of 

®vliq  *  2.546% 

Using  the  values  computed  above,  the  load-unload  response  of  the  idealized 
'saturated  porous  material  is  plotted  in  Figure  2.3  in  terms  of  total  axial 
stress,  pore  pressure,  and  effective  axial  and  radial  stresses  as  functions  of 
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volume  strain.  The  overall  behavior  during  the  uniaxial  loading  is  similar  to 
that  during  the  hydrostatic  loading,  but  with  several  significant  differences. 

•  During  unload,  liquefaction  is  assumed  to  occur  when  the  radial 
effective  stress  drops  to  zero.  Since  the  axial  effective 
stress  has  not  yet  reached  zoro,  an  instability  results.  In 
this  figure  we  have  assumed  a  strain  controlled  condition,  so 
the  instability  causes  a  sudden  drop  to  zero  in  the  axial 
effective  stress.  This  drop  in  effective  stress  results  in  a 
sudden  drop  in  total  stress  and  a  sudden  increase  in  pore 
pressure,  until  pressure  equilibrium  is  achieved  in  the 
liquefied  soil  water  mixture.  The  sudden  pore  pressure 
increase  is  due  to  the  expansion  of  the  solid  grains  caused  by 
the  sudden  drop  in  effective  stress.  The  sudden  drop  in  total 
stress  is  due  to  the  sudden  loss  of  axial  effective  stress  in 
the  material  skeleton.  The  drop  in  effective  stress  equals  the 
sum  of  the  increase  in  pore  pressure  plus  the  decrease  in 
total  stress,  according  to  the  effective  stress  law. 

•  The  actual  energy  loss  during  the  uniaxial  load-unload  cycle 
(represented  by  the  hysteresis  in  the  total  axial  stress  cycle) 
is  approximately  four  times  that  experienced  in  the  hydrostatic 
loading  to  the  seme  peak  volume  strain.  This  is  also  reflected 
by  the  higher  hysteresis  in  the  effective  axial  stress  cycle  as 
compared  to  the  effective  pressure  cycle  in  the  hydrostatic 
loading. 

•  The  pore  pressure  response  during  the  loading  and  initial  por¬ 
tion  of  the  unloading  is  identical  in  the  uniaxial  and 
hydrostatic  loadings.  At  the  onset  of  liquefaction,  however, 
the  pore  pressure  suddenly  increases  during  the  uniaxial 
unloading. 
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Effective  Volumetric  Strain,  (ev  -  j^-) 


Figure  2.1.  Idealized  simple  bilinear  skeleton  model. 
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Figure  2.3.  Undrained  uniaxial  strain  load-unload  response  of 
idealized  saturated  porous  material. 
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SECTION  3 


NUMERICAL  TREATMENT  OF  UNDRAXNEO  UNIAXIAL  LOADING 

3.1  INTRODUCTION 

As  discussed  in  the  introduction  to  Section  2,  there  are  advantages  in 
modeling  the  undrained  response  of  saturated  porous  materials  based  on  the 

properties  of  their  individual  components.  Section  2  provides  the  theoretical 

« 

basis  for  modeling  both  the  undrained  hydrostatic  and  uniaxial  response  using 
linear  approximations  of  the  pore  water  and  solid  grain  response  and  bilinear 
approximations  of  the  skeleton  response.  In  reality,  none  of  the  components 
of  saturated  porous  soils  and  rocks  are  linear,  particularly  at  the  stresses 
of  interest  in  modeling  explosive  loadings.  Basically,  all  of  the  constitu¬ 
tive  relationships  developed  in  section  2  can  be  applied  in  incremental  form 
using  numerical  techniques  to  model  the  actual  nonlinear  behavior  of  real 
earth  materials  in  undrained  loadings.  This  section  describes  the  development 
of  numerical  techniques  for  modeling  undrained  uniaxial  strain  load-unload 
behavior  of  saturated  porous  materials  and  presents  demonstration  calcula¬ 
tions.  The  numerical  model  incorporates  nonlinear  elastic  response  of  the 
solid  grains  and  pore  water  and  nonlinear  plastic  response  of  the  porous  ske¬ 
leton.  Hydrostatic  loadings,  as  discussed  in  Section  2,  can  be  treated  as  a 
special  case  of  uniaxial  loadings  where  the  constrained  moduli  equal  the  bulk 
moduli  and  Poisson's  ratio  is  0.5. 

3.2  CONSTITUTIVE  MODEL  FOR  DRAINED  POROUS  SKELETON 

The  drained  uniaxial  loading  response  of  a  porous  rock  or  soil  skeleton 
is  generally  unique  for  a  given  material  and  can  be  established  experimen¬ 
tally.  The  loading  response  can  be  digitized  and  input  as  part  of  the  skele¬ 
ton  model  for  each  particular  material.  The  unlcading  behavior  has  been  shown 
(Rlouin,  Martin  and  McIntosh,  1984)  to  be  a  function  of  the  peak  axial  stress. 
Thus,  skeleton  unloading  must  be  modeled  using  analytic  formulations  that 
describe  the  unloading  response  as  a  function  of  previous  maximum  stress. 
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Figure  3.1  shows  a  typical  uniaxial  load-unload  response  for  a  soil  ske¬ 
leton.  The  skeleton  is  strongly  hysteretic,  with  the  unload  curve  initially 
dropping  steeply,  but  breaking  back  toward  the  stress  axis  at  low  axial 
stress,  in  what  is  often  described  as  the  unloading  hook.  In  dynamic  problems 
it  is  usually  important  to  model  the  unloading  hook,  because  a  large  percen¬ 
tage  of  the  response  time  is  spent  on  this  portion  of  the  unload  curve,  and  it 
thus  has  a  strong  influence  on  late-time  displacements.  It  will  also  have  a 
strong  influence  on  the  onset  of  liquefaction  during  undrained  unloading.  The 
maximum  strain  at  the  peak  axial  stress,  a'p,  is  comprised  of  an  elastic  or 
recoverable  component,  e'e,  and  a  plastic  or  permanent  component,  e'p. 

An  expanded  view  of  the  schematic  unloading  response  is  shown  in 
Figure  3.2.  The  unloading  is  approximated  by  a  bilinear  unload  having  an  ini¬ 
tial  unload  modulus,  Mu,,  and  an  approximation  to  the  unloading  hook  of  My^. 
The  break  point  between  the  two  slopes  is  chosen  so  that  the  energy  dissipa¬ 
tion  in  the  approximation  equals  that  lost  in  the  actual  load-unload  cycle. 
This  is  achieved  by  adjusting  the  break  so  that  the  two  shaded  areas  of  Figure 
3.2  are  equal.  The  fraction,  r,  defines  the  break  point  between  the  linear 
unloading  approximations.  As  a  mathematical  convenience,  we  introduce  the 
apparent  unloading  modulus,  Nua. 

The  initial  unloading  modulus,  My,.,  can  generally  be  expressed  as  an 
exponential  function  of  the  peak  axial  stress  (Merkle  and  Oass,  1905).  It  is 
represented  by 


My-j  *  836^*4 


(3-1) 


Figure  3.3,  from  Blouin,  Martin  and  McIntosh  (1984)  shows  a  fit  to  initial 
unloading  moduli  of  Enewetak  beach  sand  given  by 


Myi  ■  347  c'p(°*65) 


(3-2) 


where  Mu^  and  o'p  are  measured  in  ksi. 
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Analysis  of  the  same  Enewetak  data  set  shows  that  the  recoverable  strain, 
e'e,  is  also  an  exponential  function  of  the  peak  stress  given  by 


(3-3) 

Figure  3.4  shows  the  fit  to  the  Enewetak  sand  data  which  is  given  by 

e'e  =  0.031  o'p(0.45)  (3-4) 

where  strain  is  expressed  in  percent  and  stress  in  psi. 

Finally,  the  break  point  between  the  two  linear  unload  approximations, 
expressed  as  the  fraction,  r,  is  approximately  constant,  as  demonstrated  in 
Figure  3,5  for  the  Enewetak  beach  sand.  An  average  r  value  of  16.2%  was 
obtained. 

From  geometric  analysis  of  Figure  3.2,  an  expression  for  the  hook  modu¬ 
lus,  MUh,  can  be  derived  which  is  given  as 

(3-5) 

where  the  apparent  unloading  modulus  is  given  by 


All  the  values  on  the  right  side  of  Equation  3-5  are  defined  in  the  previous 
equations. 

3.3  CONSTITUTIVE  MODEL  FOR  WATER 

Kim,  Blouin  and  Timian  (1986)  presented  a  composite  table  for  the 
pressure  volume  response  of  fresh  water  and  sea  water  at  pressures  from  0  to 
800  kbars.  For  calculations!  purposes  the  tangent  bulk  modulus  is  expressed 


as  a  function  of  pressure  by  fitting  the  modulus  data  at  a  temperature  of 


20°C  for  the  fresh  water  and  25°  for  the  sea  water. 

Sea  Water:  Kw  =  0,00034283tt2  +  6.1839tt  +  23806.6  (3.7a) 

Fresh  Water:  Kw  =  0.00031775tt2  +  6.2713tt  +  21746.9  (3.7b) 

for  n  <  10,000  bars 

Sea  Water:  Kw  =  9.55n  +  24758.  (3.8a) 

Fresh  Water:  Kw  =  9.49n  +  21920.  (3.8b) 

for  10,000  <  n  (  100,000  bars 

Sea  Water;  Kw  =  0.0000Q291n2  +  3.365n  +  614918.  (3.9a) 

Fresh  Water:  Kw  =  0.00000685n2  +  1.337n  +  768348.  (3.9b) 

for  n  >  100,000  bars 


where  n  and  Kw  are  expressed  in  terms  of  bars. 

Figure  3.6  parts  a  through  c,  show  the  fits  to  the  fresh  water  data. 
3.4  CONSTITUTIVE  MODEL  FOR  SOLID  GRAINS 


In  order  to  realistically  model  the  nonlinear  response  of  the  solid 
grains  to  both  the  applied  pore  pressure  and  effective  stress,  analytic 
expressions  for  the  deformation  of  solids  at  high  pressure  are  required.  High 
pressure  data  for  *any  rocks  and  minerals  show  a  linear  relationship  between 
loading  wave  velocity  and  particle  velocity  (e.g.  Allen,  1967).  Loading  wave 
velocity,  c^,  can  be  expressed  as 


cl  ■  c0  ♦  Svp 


(3-10) 


where  c0  ■  the  initial  wave  velocity  at  relatively  low 

pressure 

Vp  *  peak  particle  velocity 
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S  =  experimentally  determined  constant  relating  cl 
to  Vp  (generally  about  equal  to  1.5  for  most 
dense  rocks  and  minerals). 

Conservation  of  mass  and  momentum  on  either  side  of  the  wavefront  yield 
the  familiar  relationships  (e.g.  Chou  and  Hopkins,  1972). 


Op  =  P0CLVp 


and 


M  =  PqCL2 


(3-12) 


where  op  »  peak  axial  stress 

p0  *  initial  material  density 
H  ■  constrained  secant  modulus  defined  as 

M  - 

«p 

where  Cp  *  peak  axial  strain  corresponding  to  the  peak 

stress  Op 

Substitution  of  Equation  3-10  into  3-11  gives 


(3-13) 


Op  ■  P0C0Vp  ♦  PoSVp2 


(3-14) 


and  solving  for  peak  particle  velocity  as  a  function  of  peak  stress  yields 


where 


(3-16) 
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(3-16) 


f(ffp)  =  (P02c02  +  4P0Sap)*  -  P0C0 


Substitution  of  Equation  3-15  into  3-12  gives 


M  =  F(<rp) 


where 


(3-17) 


F(ap)  =  p0c02  +  c0f(o) 


(3-18) 


The  constrained  tangent  modulus,  Mt,  used  in  the  numerical  model  is 
defined  as  the  slope  of  the  stress  strain  curve  by 


N*  ■ 


dg 

de 


Equating  3-17  and  3-13  and  solving  for  €p  yields 


(3-19) 


(3-20) 


Differentiating  Equation  3-20  with  respect  to  op  and  inverting  gives  the 
constrained  tangent  modulus  as 


M* 


F(ffp)  “  ffpF'(Pp) 


(3-21) 


Differentiating  Equations  3-18  and  3-16  with  respect  to  <rp  yields 
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(3-22) 


F'  (ffp)  =  Cgf  '  (Cfp) 


+ 


f(op)f'(cp) 

2?0 


and 


f'(ap) 


2P0S _ 

"  (P02c o2  +  4p0Sffp)< 


(3-23) 


Equations  3-15  through  3-23  can  be  used  to  define  the  high  pressure 
constrained  stress  strain  and  modulus  relationships  for  the  solid  grains.  A 
numerical  example  is  shown  in  Figure  3.7  using  the  following  material 
parameters : 


p0  *  5.14  lb-S2/ft*  (2650  kg/«3) 
c0  -  16735  ft/s  (5100  a/s) 

S  -  1.5 

Substituting  c0  into  Equation  3-12  gives  an  initial  secant  modulus  of 

M  ■  10  x  106  psi  (690  kbars) 

As  explained  in  Section  2,  the  volumetric  relationships  for  the  solid 
grains  should  be  specified  in  terms  of  the  bulk  modulus,  Kg,  rather  than  in 
terms  of  the  constrained  modulus.  At  high  pressures,  the  shear  strength  of 
the  grain  materials  becomes  insignificant  compared  to  the  applied  stress  and 
the  materials  tend  to  behave  like  fluids.  At  these  pressures,  the  bulk 
tangent  modulus  equals  the  constrained  tangent  modulus  with  Poisson's  ratio 
equal  to  0.5.  Beneath  some  threshold  stress,  Pp,  Poisson's  ratio  begins  to 
decrease  from  0.5  at  Pp  to  an  initial  value  of  Poisson's  ratio,  w0,  at  a  low 
value  of  mean  stress.  Me  have  used  the  simple  relationship  depicted  in  Figure 
3.8  to  approximate  the  influence  of  mean  stress  on  Poisson’s  ratio  and  the 
ratio  of  bulk  modulus  to  constrained  modulus,  g(p)  is  related  as 
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(3-24) 


Kg  =  g(p)Mt 


where  g(p)  is  defined  as 


for  p  <  pb. 
For  p  > 


9(P) 


2  (1  -  2u0)  fi_  *  (1  +  *>o) 

3  (1  -  i>o>'  Pb  3TT  -155) 


fl(p)  *  1 


(3-25) 


(3-26) 


Poisson's  ratio  can  be  computed  as  a  function  of  the  Modulus  ratio,  g(p),  at  a 
given  pressure  as 


v  * 


1  ♦  3g(p) 


(3-27) 


Figure  3.6  shows  both  the  Modulus  ratio  and  Poisson's  ratio  plotted  as  a  func¬ 
tion  of  pressure  for  typical  solid  grains  having  an  initial  Poisson's  ratio  of 
0.2  and  a  threshold  pressure,  pb,  of  5  kbars. 


3.5  MODELING  PRELOAD  CONDITIONS 


Dynamic  undrained  loadings  of  in  situ  saturated  Materials  do  not  start 
froM  a  condition  of  zero  stress  because,  except  right  near  the  surface,  the 
materials  are  under  some  initial  in  situ  stress  conditions.  The  effective 
vertical  stress  will  equal  the  stress  ieposed  by  the  buoyant  weight  of  the 
overburden,  though  this  nay  be  modified  by  local  tectonic  stress  conditions. 
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The  pore  pressure  will  normally  be  a  function  of  the  water  table  depth,  which 
will  also  affect  the  buoyant  overburden  load  and  effective  stresses.  Thus, 
initial  in  situ  effective  stress  and  pore  pressure  conditions  are  infinitely 
variable  relative  to  one  another  and  must  be  treated  as  mutually  independent. 
In  order  to  compute  the  undrained  response  to  uniaxial  strain  loadings  of  in 
situ  material,  a  methodology  has  been  developed  to  account  for  the  influence 
of  the  in  situ  stresses  and  pore  pressures  on  the  in  situ  material  response. 


It  is  helpful  to  use  a  laboratory  test  procedure  analogy  to  develop  the 
governing  in  situ  stress-strain  relationships.  We  can  establish  three  sample 
states: 


1.  The  virgin  stress-free  (unstrained)  state  prior  to  application 
of  any  effective  stress  or  pore  pressure;  we  assume  a  total  volume 
of  the  saturated  sample,  V0,  of  unity; 

V0  -  1  (3-28) 

The  corresponding  porosity,  h0,  is  defined  as 

fi0  «  y.q=~..vgq  .1-^22  (3-29) 

v0  v0 

where  Vgo  *  volume  of  solid  grains  in  the  virgin  unstrained  state. 

Equation  3-29  can  be  rewritten  to  give  the  grain  fraction  in  terms 
of  porosity  as 


yao  =  1  -  n 
*o 


o 


(3-30) 


2.  The  in  situ  stress  state,  where  the  in  situ  effective  stresses  and 
pore  pressures  are  either  known  or  can  be  defined.  The  in  situ 
stresses  and  pore  pressures  are  governed  by  the  depth  of  overburden 
and  depth  of  the  water  table.  The  notation  for  the  in  situ  state  is 
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»  in  situ  pore  pressure 

<r'j  =  in  situ  effective  vertical  (axial)  stress 
Vw^  =  volume  of  pore  water  in  the  in  situ  stress 
condition 

3.  A  set  of  hypothetical  parameters  which  preserve  the  in  situ  material 
mass  in  the  zero  stress  condition.  In  addition  to  Vg0,  we  have 

VWo  =  volume  of  pore  water  in  the  in  situ  material, 
but  at  zero  pore  pressure 

Define  the  total  volume  of  in  situ  material  in  the  zero  stress 
condition,  V0»  as 


Vo  =  Vgo  +  Vwo  (3-31) 

Define  the  ratio  of  Vwo  to  V0  as  a  fictitious  porosity,  n0, 
according  to 

n0  ■  ^  (3-32) 

V0 

The  parameters  in  this  third  state  are  defined  for  mathematical 
convenience.  Since  the  pore  pressure  and  the  effective  stress  in  the 
in  situ  condition  are  specified  independently,  the  volumes  of 
in  situ  material  when  taken  to  the  zero  stress  state  do  not  gen¬ 
erally  constitute  a  compatible  saturated  material.  The  volume 
of  zero  stress  pore  water  will  generally  be  too  small  or  too  large 
to  fill  the  available  pore  space. 

Equation  2-15  can  be  rewritten  in  incremental  form  as 

dE  -  (3-33) 

where  the  modulus  values  represent  the  tangent  moduli  at  the  stress  level  of 
interest  and  the  total  strain  increment  de  is  in  terms  of  the  virgin  state. 
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The  final  term  dn/3Kg  represents  the  axial  component  of  strain  created  by  the 
pore  water  compression  of  the  solid  grains.  It  can  be  expressed  as 

d®gir  =  jij^  (3-34) 

where  egn  is  the  volumetric  strain  of  the  solid  grains  due  to  the  pore 
water  pressure. 

The  term  3Ks/Hs  represents  the  apparent  softening  of  the  skeleton  due  to  the 
fact  that  the  uniform  decrease  in  grain  volume  generated  by  the  pore  pressure 
must  be  redistributed  as  an  axial  volume  decrease  in  the  uniaxial  loading. 

The  term  do'/Hs  is  the  component  of  strain  due  to  the  effective  stress  and 
can  be  expressed  as 


d2'  -  —  (3-35) 

"s 

Using  Equations  3-34,  3-35  and  2-25,  Equation  3-33  can  be  rewritten  in  terms 
of  incremental  strain  as 


« ■ «' +  iltVi,  <fc8*  .  <3-36> 

To  obtain  the  total  strain  at  the  in  situ  stress  state,  we  integrate 
Equation  3-36  to  give 


ei  »  «i'  +  3(1  !  l\  sgirt 


(3-37) 


where  the  terms  C('  and  represent  integrations  of  the  nonlinear  strain 
relationships  with  respect  to  stress. 

The  total  strain  in  the  in  situ  condition  is  defined  as 

ii  «  1  -  p-  (3-38) 
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and  the  ratio  of  the  in  situ  total  volume  to  the  initial  volume  in  the  virgin 
state  is  obtained  by  combining  Equations  3-37  and  3-38  to  give 


-U.+.Hl 
3(1  -  v) 


€grri 


(3-39) 


The  total  strain  in  the  solid  grains  in  the  in  situ  condition,  cgi,  is 
defined  as 


egi  «  1  -  M  (3-40) 

VgO 

where  Vgi  =  the  in  situ  grain  volume 

This  strain  is  a  function  of  the  pressure  on  the  grains,  pg^,  which  is  given 
by 


(3-41) 


where  n^  ■  in  situ  porosity 

p'l  •  effective  mean  stress  in  the  porous  skeleton 
p'i/(l~n^)  «  pressure  on  the  solid  grains  due  to  effective 
stresses  in  the  skeleton. 


The  effective  mean  stress  in  the  skeleton  can  be  determined  from 


P‘i  8  1/3  (d’i  ♦  20’,., i 


where  a'r^  •  the  in  situ  radial  effective  stress 

From  equations  2-22,  2-23  and  3-42  we  obtain 


P'i 


-  M  ♦-») 

3(1  -  u) 


*’i 


2(.l.--M  *s 
3(1  -  y)  K£ 


«i 


(3-42) 


(3-43) 
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where  the  bulk  moduli  values  are  the  secant  moduli  at  the  in  situ  stress 
level. 

The  total  strain  in  the  grains  in  the  in  situ  condition  can  now  be 
obtained  by  using  the  total  pressure  on  the  grains  from  Equation  3-41  to  com¬ 
pute  the  strain  from  the  grain  properties  described  in  Section  3.4.  The 
volume  of  the  grains  in  the  in  situ  state  relative  to  the  volume  in  the  virgin 
state  is  obtained  by  rewriting  Equation  3-40, 

(3-44) 

Finally,  the  volume  strain  in  the  pore  water  in  the  in  situ  state  is 
given  as 

e**  *  1  -  M  (3-46) 

vwo 

and  can  be  directly  obtained  from  the  pressure  volume  re’ationships  for  water 
described  in  subsection  3.3,  Equation  3.45  is  rewritten  as 

(3-46) 


The  total  compressed  volume  of  material  in  the  in  situ  state  is  given  by 

(3-4?) 

The  calculational  procedure  used  to  compute  the  initial  in  situ  state  is 
outlined  below. 

1.  let  the  virgin  volume,  V0,  equal  unity  using  Equation  3-28. 


33 


2.  Compute  the  solid  grain  volume,  Vgo,  from  Equation  3-30. 

3.  Compute  the  solid  grain  volume  under  the  in  situ  stress  conditions, 

Vgi,  from  Equation  3-44. 

4.  Compute  the  total  in  situ  volume  from  Equation  3-39. 

5.  Compute  the  in  situ  pore  water  volume  from  Equation  3-47. 

6.  Compute  the  in  situ  pore  water  volume  at  zero  pressure,  Vwo* 

from  Equation  3-46. 

7.  Compute  the  total  volume  of  in  situ  material  in  the  2ero  stress 
state,  V0,  from  Equation  3-31. 

8.  Compute  the  fictitious  porosity,  n0,  from  Equation  3-32. 

9.  Compute  the  in  situ  porosity,  n^,  from 


10.  Compute  the  in  situ  volume  strain,  c^,  from  Equation  3-36. 

3.6  NUMERICAL  MODELING  OF  UNIAXIAL  UNDRAINED  RESPONSE 
3.6.1  Computational  Algorithm  NKOCP 

In  order  to  realistically  predict  the  response  of  saturated  porous 
materials  to  undrained  loadings,  the  nonlinear  response  characteristics  of  the 
individual  components  must  be  accurately  modeled.  This  requires  an  incremen¬ 
tal  numerical  technique  which  can  incorporate  all  important  behavioral  charac¬ 
teristics  of  the  solid  grains,  pore  water  and  most  importantly  the  material 
skeleton.  The  algorithm  NKOCP  fully  describes  the  load-unload  response  of 
fully  saturated  porous  materials  in  uniaxial  strain  conditions.  NKOCP  uses 
the  fully  coupled  material  model  described  by  Kim,  Blouin  and  Tinian  (1986). 

It  can  also  be  used  to  simulate  hydrostatic  loadings  by  substituting  tha  bulk 
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moduli  for  the  constrained  moduli,  and  setting  Poisson's  ratio  to  0.5.  NKOCP 
is  an  advanced  version  of  KOCP,  described  by  Kim,  Blouin  and  Timian  (1986) 
with  the  added  capabilities  of: 

•  modeling  the  response  of  materials  starting  from  an  arbitrary 
specified  in  situ  stress  condition; 

•  modeling  the  unloading  response  including  liquefaction. 

This  improved  algorithm  incorporates  the  unloading  model  for  the  skeleton 
described  in  subsection  3.2,  the  solid  grain  model  described  in  3.4  and  the 
preload  capabilities  described  in  3.5.  Details  of  the  computational  proce¬ 
dures  are  described  by  Kim,  Blouin,  and  Timian  (1986),  Section  5. 

3.6.2  NKOCP  User's  Manual 

NKOCP,  listed  in  Appendix  A,  uses  input  from  Tapes  55  and  8.  Tape  55 
includes: 

•  the  digitized  skeleton  stress-strain  loading  response  (usually 
obtained  from  drained  laboratory  tests), 

•  parameters  describing  the  drained  skeleton  unloading  charac¬ 
teristics, 

•  in  situ  parameters 

•  solid  grain  properties 

•  general  input/output  control  parameters 

Tape  8  contains  a  table  of  K0  vs.  effective  vertical  strain  pairs.  The 
compressibility  of  the  pore  water  from  subsection  3.3  is  built  in.  All  inputs 
are  in  free  format. 

DESCRIPTION  OF  INPUT  FILES 

TAPE  55 
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Card  1  -  General  Parameters 


NF,  NW,  EMAX,  NDIV,  NDIVU,  NDIVM,  LPRINT 


NF  =  0  pressure  in  psi 
=  1  pressure  in  bars 

=  2  pressure  in  psi  with  constant  water  modulus,  Kw,  of 
0.3  x  10°  psi  and  constant  solid  grain  modulus.  Kg,  of 
5  x  106  psi . 

NW  =  0  fresh  water 
»  1  sea  water 

EMAX  =  Maximum  total  strain  at  the  start  of  unload 


NDIV  «  Number  of  interpolation  divisions  between  each  consectutive 
data  point  on  the  loading  stress-strain  curve 

NDIVU  *  Number  of  equal  strain  decrements  between  e'p  and  EMAX 
used  to  compute  unloading  prior  to  liquefaction 

NDIVM  b  Number  of  equal  strain  decrements  between  0  and  the  point 
of  liquefaction  used  to  compute  unloading  following  lique¬ 
faction 


LPRINT  =  Output  printing  interval,  e.g.  for  LPRINT  »  10  every  tenth 
computed  point  will  be  output. 


Card  2  -  In  Situ  Characteristics 


PF,  SE,  V,  POR 


PF  c.  in  situ  pore  water  pressure 
SE  «  In  situ  effective  vertical  stress 


V  a  Poisson's  ratio  of  in  situ  skeleton 


36 


POR  =  Virgin  porosity  prior  to  any  loading. 
Card  3  -  Skeleton  Unloading  Characteristics 
Al,  A2,  A3,  A4,  R,  XKOI,  XKOH 


Al,  A2  =  Recoverable  unloading  strain  parameters  from  Equation  3-3 

A3,  A4  =  Parameters  for  initial  unloading  constrained  modulus  from 

Equation  3-1 

R  =  Stress  fraction  at  the  break  point  of  the  two  unloading 
slopes  described  in  Figure  3.2 


XKOI  =  Value  of  K0  during  the  initial  unloading 

XKOH  =  Value  of  K0  during  unload  along  the  second  (hook)  portion 

of  the  bilinear  unloading  slope 

Card  4  -  Solid  Grain  Properties 


RO,  CO,  GVO,  SS,  PB 


RO  e  initial  grain  density,  lb-s2/in4 

co  «  Initial  wavespeed,  in/s 

GVO  «  Poisson's  ratio  at  low  stress  from  Figure  3.8 

SS  «  Slope  of  propagation  velocity  to  particle  velocity 
relationship,  S 

PB  s  Threshold  pressure,  psi,  where  the  tangent  constrained 
modulus  becomes  equal  to  the  tangent  bulk  modulus 
from  Figure  3.8 
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Card  5 


Digitized  Virgin  Skeleton  Loading  Stress-Strain  Response 


NPOINT 


e’l,  o’! 
e'2.  o' 2 


£'n»  flr'n 


TAPE  8 


NPOINT  cards  with  effective  vertical  strain 
and  stress  pairs 


Card  1  -  Kq  Values 


for  Skeleton  During  Loading 


NKO 


s'l.  ^1 

2>  ko2 


c '  n »  *on 


NKO  cards  with  effective  vertical  strain  and  K0 
pairs,  e.g.  between  c'j  and  e*2  the  value  Ko2 
is  used 


DESCRIPTION  OF  OUTPUT  FILES 

Output  is  on  three  tapes,  6,  3,  and  66.  Tapes  6  and  3  contain  load-unload 
response  prior  to  liquefaction.  Tape  66  contains  the  response  following 
liquefaction. 
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w  v»y»>-«yv-  jmw  3^* 


TAPE  6 


e,  o\,  it,  a,  LF 

I  I  II  l 


e  =  Total  volume  strain  (%) 
a'v  =  Effective  vertical  stress 
jr  s  Pore  water  pressure 
a  =  Total  vertical  stress 

LF  a  o  indicates  stress  state  on  virgin  loading  curve 

*  1  indicates  stress  state  on  the  initial  unloading  curve 
-  2  indicates  stress  state  on  the  hook  portion  of  the 
unloading  curve 

s  4  indicates  stress  state  post  liquefaction,  mixture 
response 

TAPE  3 

e,  o*h,  Mf*  n 

III  I  I 

I  I  I  I  I 


e  *  Total  volume  strain  ( % ) 
o*h  ■  Effective  horizontal  stress 
0^  *  Total  horizontal  stress 
Mf  *  Fully  coupled  constrained  tangent  aodulus 
n  *  Current  porosity 
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TAPE  66 


€,  ft,  Km,  n 


e  »  Total  volume  strain  (X) 
it  =  Pore  water  pressure 
Km  3  Undrained  tangent  mixture  modulus 
n  »  Current  porosity 

3.6.3  Example  Problem 

An  example  problem  Is  documented  In  Appendix  B.  This  example  problem 
duplicates  the  sample  problem  described  In  Section  2.4.  Appendix  B  Includes 
Table  B.l  for  the  listings  of  Input  files,  Table  B.2  for  the  listings  of  out¬ 
put  files,  and  Figure  B.l  for  the  load-unload  response  of  the  example  problem. 


40 


Figure  3.1.  Schematic  skeleton  stress-strain  behavior  during  drained  uniaxial 
strain  loading. 


Figure  3.2. 


z 1 


Variable  bilinear  approximation  of 
constrained  unloading. 
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Figure  3.4.  Recoverable  strain  as  a  function  of  peak  stress  for  Enewetak  beach  sand. 
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Figure  3.8.  Poisson's  ratio  and  modulus  ratio  g{p)  as  a  function 
of  pressure  for  solid  grains. 
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SECTION  4 

THEORETICAL  ANALYSIS  OF  WAVE  PROPAGATION 
IN  SATURATED  POROUS  MEDIA 


4.1  INTRODUCTION 

The  fully  coupled  field  equations  presented  by  Kim,  BLouin  and  Timian, 
1986,  Section  7,  were  used  to  derive  an  analytical  closed-form  solution  for 
wave  propagation  velocity  in  fully  saturated  porous  media.  This  solution  is 
highly  versatile  in  that  it  considers  compression  of  the  solid  grains, 
compression  of  the  pore  water,  deformation  of  the  porous  skeleton,  and  spatial 
damping  and  can  be  used  to  compute  wavespeeds  of  the  first  and  second  kind  and 
damping  coefficients  in  various  geologic  materials.  This  solution  provides  a 
means  of  verifying  the  multiphase  code  calculations  and  also  provides  a  means 
of  analyzing  the  influence  of  material  property  variations  on  wavespeed  and 
attenuation.  Subsections  4.2  and  4.3  describe  the  development  of  the  closed 
form  solution  for  the  wavespeeds  of  the  first  and  second  kind. 

4.2  GENERAL  RELATIONSHIPS 

The  behavior  of  saturated  porous  media  can  be  described  by  two 
equations  of  motion  and  two  supplementary  equations.  The  first  supplementary 
equation  is  that  for  continuity  of  flow  and  the  second  describes  the  effective 
stress  law.  The  first  equation  of  motion  prescribes  the  motion  of  the  bulk 
mixture  and  the  second  defines  the  motion  of  the  pore  fluid  relative  to  the 
porous  skeleton. 

4.2.1  Governing  Equation  for  Bulk  Mixture 

The  differential  equation  governing  the  bulk  mixture  is  given  by 

*ij,j  *  U  -  «)  Pg  u<  ♦  n  pf  U<  (4-1) 

°ij,j  t0**l  stress  gradient  applied  to  an  infinitesimal  element  of 

saturated  material  at  some  given  time,  Ofj#j  is  expressed  in  tensor  notation 


51 


and  represents  the  stress  gradient  in  each  of  3  mutually  perpendicular  coor¬ 
dinates  (e.g.  see  Mendelson,  1968).  For  instance,  in  the  x  direction, 


J .  j  “  a 


*^xx  yy  Bcr 


xz 


x 


=  (1  -  n)  Pg  Uy  +  n  p^  Uy 


(4-2) 


The  term  (1  -  n)pg  is  the  mass  of  the  soil  skeleton  per  unit  volume  of 
saturated  material,  where  n  is  the  porosity  and  pg  is  the  mass  density  of  the 
solid  grains,  u-j  is  the  dispalcement  of  the  skeleton  in  the  i  direction  and 

u-j  is  the  acceleration  of  the  skeleton  in  the  i  direction.  The  term  n  p^  is 
the  mass  of  the  pore  fluid  per  unit  volume  of  saturated  material  where  p^  is 
the  mass  density  of  the  pore  fluid.  U-j  is  the  absolute  displacement  of  the 

pore  fluid  in  the  i  direction  and  U-j  is  the  absolute  acceleration  of  the  pore 
fluid  in  the  i  direction. 


The  bulk  mass  density  of  the  saturated  material,  p,  is  given  by 

P  =  (1  -  n)pg  +  n  pf  (4-3) 

Substitution  of  the  value  for  (1  -  n)pg  from  Equation  4-3  into  Equation  4-1 
gives 


<Mj,j  ■  (P  -  n  pf)u1-  ♦  n  pf  (4-4) 

A  term  w j  is  introduced  which  is  the  apparent  fluid  displacement  in  the  i 
direction  relative  to  the  soil  skeleton  and  is  given  by 

w-j  =  n(U j  -  Uj)  (4-5) 

In  seepage  problems,  w-j  is  referred  to  as  the  discharge  displacement.  It 
describes  the  discharge  of  fluid  through  a  soil  mass  of  unit  area.  The 

discharge  velocity,  or  apparent  relative  velocity,  w-j,  between  the  soil  par¬ 
ticles  and  pore  water  is  the  velocity  of  water  in  a  discharge  duct  of  unit  area 
needed  to  maintain  the  actual  relative  velocity  in  the  porous  soil  of  the  same 
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unit  area.  The  actual  relative  velocity  between  the  skeleton  and  the  pore 

water  is  given  by  w-j/n.  Finally,  w,  is  the  apparent  relative  acceleration 
between  the  soil  skeleton  and  pore  water  given  by 

iii  =  n(Ui  -  u,-)  (4-6) 


Equation  4-6  is  an  approximation  of  the  actual  relative  acceleration  because 
the  porosity  is  not  constant,  as  assumed  here,  but  is  in  reality  changing  with 
time.  This  is  just  one  of  several  simplifying  assumptions  made  in  the  linear 
analysis  which  make  this  theoretical  development  tractable.  Equation  4-4  can 
be  expressed  in  terms  of  the  apparent  relative  fluid  acceleration  as  simply 


aij,j  s  Pui  *  Pf*i 


(4-7) 


4.2.2  Governing  Equation  for  Pore  Fluid 

The  equation  of  motion  governing  the  pore  fluid  is  derived  from 
Biot's  work  (1956).  Biot  expressed  the  pore  pressure  gradient,  n, (,  as 


n,i  «  pf  U*  ♦  Di  (4-8) 

where  the  pore  pressure  gradient  is  expressed  in  tensor  notation. 

For  example,  the  gradient  in  the  x  direction  is  given  by 

x,x  *  Pf  ^  ^x  (4-9) 

The  term  pfU*  is  the  inertial  force  per  unit  volume  of  pore  fluid.  repre¬ 
sents  the  viscous  friction  force  between  the  pore  fluid  and  the  soil  skeleton 

per  unit  volume  of  pore  fluid.  Solving  equation  4-6  for  and  substituting 
into  Equation  4-8  gives 


*'i  “  “n  "i  ♦  Pf  «i  ♦  °i 


(4-10) 
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For  dynamic  laminar  flow,  Biot  (1956)  derived  the  exact  expression  for  the 
viscous  friction  term,  Df,  for  both  circular  and  flat  ducts  which  is  given  by 

□i  =  Tf  ^  Wi  (4-11) 

where  yf  is  the  unit  weight  of  pore  fluid  and  F(jc)  is  the  viscous  friction 
correction  factor,  and  k  is  the  coefficient  of  permeability.  F(k)  is  a 
complex  function  given  by 


F(k)  =  fiOc)  +  i  f20c) 


(4-12) 


The  magnitudes  of  the  real  part,  fj(ic),  and  imaginary  part,  f2(ic)  are  plotted 
in  Figure  4.1  as  a  function  of  the  nondimensional  parameter,  jc,  which  is 
defined  as 


k  *  r  kq  (4-13) 

/ 

where  ic©  ■  ((~)  w  k|*  (4-14) 

In  Equation  4-14,  g  is  the  gravitational  accleeration,  w  is  the  excitation 
frequency,  and  the  factor,  r,  is  the  constant  which  is  dependent  on  the  shape 
of  flow  path.  For  Generalized  Darcy's  Flow  (also  called  Poiseuille  Flow), 
fj(x)  and  f2U)  are  independent  of  tc  and  r  »  0.  For  the  circular  duct,  r 
is  unity.  And  for  the  flat  duct,  the  value  of  r  is  approximately  equal  to 


For  the  purpose  of  simplifying  expressions,  we  define  a  modified  coef¬ 
ficient  of  permeability  as 


ic 


(4-15) 


Substituting  Equation  4-15  into  4-11  gives 


D 


i 


a 


Vf 

k 


W-i 


(4-16) 
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Note  that  as  shown  in  Equation  4-12,  that  k  is  a  complex  function  of  jc. 

An  approximate  form  of  Equation  4-16  was  developed  by  Kim  and  Blouin 
(1984)  and  is  given  by 


Di  °  IT  +  YT  r  "i  (4-17) 

where  Yf  =  unit  weight  of  pore  fluid 

r  =  empirical  mass  increment  factor 

The  theoretical  value  of  r  is  1/3  for  the  circular  duct  and  1/5  for  the  flat 
duct.  Equation  4-17  is  a  good  approximation  when  the  nondimensional  parameter, 
k,  is  less  than  2. 

For  mathematical  convenience,  and  to  allow  us  to  compare  the  exact  and 
approximate  solutions,  both  solutions  are  incorporated  in  a  single  expression 
for  the  fluid  friction  given  by 

■  If  r  W{  (4-18) 

k  n 

In  using  Equation  4-18,  it  should  be  noted  that  for  the  approximate  solution 
r  *  0  and  k  =  k,  and  for  the  exact  solution  r  ■  0  and  it  *  it.  Substitution  of 
Equation  4-18  into  Equation  4-10  gives 


rr,i  »  (1  ♦  r)  +  pfii*  ♦ 


(4-19) 


4.2.3  Continuity  Equation  of  Flow 

The  continuity  equation  for  pore  fluid  flow  is  derived  from  mass  con¬ 
servation  relationships.  The  volumetric  strain  of  the  pore  fluid,  ef,  is 
given  by 


dpf 

def  ■  — “  »  Cf  dir 


(4-20) 
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where 


Cf  =  pore  fluid  compressibility 
n  =  pore  fluid  pressure 

The  volume  strain  of  the  solid  grains,  eg,  is  given  by 

de9  =  "  =  c9  drr  +  1?^  dp'  {4'21) 

where 


Cg  =  bulk  compressibility  of  solid  grains 
p'  =  effective  mean  pressure 

The  dry  density,  pj,  is  given  by 


Pd  =  5?  =  (1  "  n)pg  (4“22) 

where  mg  is  the  mass  of  the  solid  grains  in  skeleton  volume  V*.  The  change  in 
dry  density  is  given  by 

dpd  =  -Pd  dev  (4-23) 

where  ev  is  the  volumetric  strain  of  the  skeleton.  Differentiating  Equation 
4-22  with  respect  to  n  and  pg  gives 

dpd  ■  (1  -  n)  dpg  -  pgdn  (4-24) 

Equating  4-23  and  4-24  yields 


dc  a 

v  ‘  l-n  Pg 


Conservation  of  mass  for  the  pore  fluid  within  a  specified  initial 
volume  of  saturated  porous  material  is  given  by 

n  pf  Vt  =  n'pf'Vt'  (4-26) 
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where,  as  illustrated  in  Figure  4.2,  the  terms  to  the  left  of  the  equal  sign 
represent  the  fluid  mass  under  the  initial  conditions  and  the  terms  to  the 
right  represent  the  same  fluid  mass  under  deformed  conditions.  Equation  4-26 
may  be  expressed  in  infinitesimal  incremental  form  as 

npfVt  =  (n  +  dn)  (pf  +  dpf)  (1  +  dep)  Vt  (4-27) 

where 

€p  =  volumetric  diffusion  of  pore  fluid  as  depicted  in  Figure 
4.2. 


Solving  Equation  4-27  for  dep  and  discarding  second  order  terms  yields 


(4-28) 


Equation  4-28  is  combined  with  Equation  4-25  to  yield 
(1  -  n)d€v  ♦  ndep  +  (1  -  n)  ♦  n  —  »  0 

Pg  Pf 

Combining  Equations  4-20  and  4-21  with  4-29  gives 

n(d£p  -  dev)  ♦  dev  -  ^  dn  -  cfldp'  «  0 


(4-29) 


(4-30) 


where  Ka  is  the  bulk  modulus  of  the  solid/fluid  mixture  which  is  expressed  by 


*m  B 


ncf  +  (l  -  n)Cg 


(4-31) 


The  change  in  effective  mean  pressure  is  given  by 

dp'  »  Ks  (dcv  -  Cgdx)  (4-32) 

Substituting  Equation  4-32  into  4-30  gives 
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(4-33) 


n(d€p  -  dev )  +  (1  -  CgK§)  dey  +  (Cg^Kg  - 


J-)  dn  =  0 


4.2.4  Effective  Stress  Law 


From  Terzaghi's  effective  stress  law,  the  change  in  total  stress  is 
expressed  as  the  sum  of  the  change  in  effective  stress  and  the  change  in  pore 
water  pressure  as 


da.jj  =  do'-jj  +  6-f j  dn 


where 

a\ j  =  total  stress 
o'-jj  *  effective  stress 
»  Kronecker's  delta 
fiij  »  0  if  i  i»  j 
fiij  ■  1  if  i  »  j 


(4-34) 


4.3  PROPAGATION  VELOCITY  OF  PLANE  COMPRESSIONAL  WAVES  IN  SATURATED  ELASTIC 
POROUS  MEDIA 


Using  the  four  general  relationships  derived  in  Section  4.2, 
equations  for  dynamic  response  to  uniaxial  strain  loadings  can  be  easily 
determined.  These  have  broad  general  application  to  explosive  loadings  which 
typically  generate  planar  or  spherical  compressions!  waves  in  saturated  porous 
media. 


4.3.1  Governing  Equations  for  Uniaxial  Strain  Loadings 

For  plane  wave  propagation  in  the  vertical  direction  Equation  4-7  for 
the  bulk  mixture  can  be  expressed  as 
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(4-35) 


Likewise,  Equation  4-19  governing  motion  of  the  pore  fluid  is  expressed  as 


pf  a? 


du2  + 


*,  +  —  (! 
2  n 


dt 


(4-36) 


The  continuity  equation  of  flow,  Equation  4-33,  can  be  differentiated 
with  respect  to  time  to  give 


n(eF  -  cv)  +  (1  -  cgKs)ev  +  (CqZKj.  -  j^)*  =  0  (4-37) 

From  Figure  4.2  it  can  be  seen  that  the  strain  difference,  eF  -  ev,  represents 
the  extent  of  the  pore  water  diffusion  relative  to  the  solid  skeleton.  Thus, 
the  apparent  relative  velocity  of  the  pore  water  can  be  expressed  as 

♦ 

d  W 

n(£F  -  *v>  *  a f  (4“38) 

The  volumetric  strain  rate  of  the  porous  skeleton,  £v,  can  be  expressed  as 

8  dUv 

-  3t  ST  <4‘39> 

L  duJ 
a  3z  8t 

“  dZ 

Substituting  Equations  4-38  and  4-39  into  4-37  gives 
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(4-40) 


8Wy  9UZ  •  1  Ajr 

a T  +  (1  '  cgKs)  aT  +  cg2|<s  "  at  =  0 


The  effective  stress  law  from  Equation  4-34  can  be  expressed  for 
uniaxial  strain  as 


dcz  =  d<x2'  +  dw  (4-41) 

The  effective  vertical  stress  is  a  function  of  the  vertical  strain,  e2,  and 
the  pore  fluid  pressure,  tt,  applied  to  the  solid  grains  according  to 

do2'  =  Mg  (de2  —  CgKgdTT)  (4—42) 

where  Ms  is  the  constrained  modulus  of  the  porous  skeleton  (see  Blouin  and 
Kim,  1984).  Substitution  of  Equation  4-42  into  4-41  yields 

da2  *  Msdez  +  (1  -  CgKg)dfT  (4—43) 

Differentiating  Equation  4-43  with  respect  to  time  gives 

&z  *  Mg£z  ^  (1  “  CgKg)A  (4-44) 

For  uniaxial  strain 

*  £v  (4-45) 

Combining  Equations  4-39,  4-44,  and  4-45  gives 


(1 


an  aoz 
coKs>  at  "  aT 


o 


(4-46) 
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4.3.2  Wave  Velocity  for  Harmonic  Excitation 

For  vertically  propagating,  one  dimensional  oscillatory 
(sinusoidal)  motion  with  the  pore  fluid  and  porous  skeleton  in  phase.  The 
response  of  the  saturated  porous  medium  may  be  expressed  as 


uz  =  uzo  eaz  e^w(*  *  d) 

(4-47) 

wz  =  wzo  ea7‘  eiw(t  ♦  §) 

(4-48) 

n  =  eaz  e*w(*  ♦  c) 

(4-49) 

crz  =  czo  ®az  e*w(*  *  s) 

(4-50) 

where  c  is  the  wave  velocity  and  a  is  the  damping  constant  for  a 

tion  frequency  w.  For  mathematical  simplicity  we  define  c*  such 

given  excita- 

that 

*  .  i  -  iS 

c*  c  Ci> 

(4-51) 

Substitution  of  Equation  4-51  into  4-47  through  4-50  yields 

uz  »  uzo  e*w(*  *  §*) 

(4-52) 

wz  ■  wzo  e*w(*  4  l‘ ) 

(4-53) 

n  >-  n0  e*w(t  *  I1  * 

(4-54) 

Oz  m  <?zo  *  <5*  V 


(4-56) 


Differentiating  Equations  4*52  through  4-55  with  respect  to  time  yields 


auz 

aT  =  1  w  uz 


(4-56) 


3*, 

aF  =  1  w  w* 


(4-57) 


an 

at 


=  i  w  n 


(4-58) 


3oz 

a r  a  s  w  az 


(4-59) 


Differentiating  Equations  4-52  through  4-55  with  respect  to  z  yields 


^z  iu  . 

sT ’?4i 


(4-60) 


<h&Z  i«  • 
aT ’F"* 


(4-61) 


an  iu 
dz  c' 


(4-62) 


dz 


(4-63) 


Substituting  the  differential  Equations  4-56  through  4-63  into  the  governing 
equations  for  the  uniaxial  strain  loadings,  4-35,  4-36,  4-40  and  4-4$  yields 
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Equation  4-64  can 

to  zero.  A  modest  amount 

be  solved  by  setting  the  determinant  of  the  matrix 

of  algebraic  manipulations  yields 

a4(c*)4  -  b2(c')2  +1=0 

(4-65) 

where 

*4 . 

(4-66) 

b2  _  -b3  ♦  »«i 

(4-67) 

and 

bj  B  bj  -  b2f2(»0 

(4-68) 

b2  =  b2f  j(ic) 

\4— 69 ) 

b3  =  B3  +  b4f2(x) 

(4-70) 

b4  =  b4f  ^  (fc) 

(4-71) 

and 

*>1 

3  P  *2  *3  ”  Pf 2  a3 

(4-72) 

b2 

a  Ml  83 

(4-73) 
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b3  =  2  Pf  a4  -  s>2  a42  ~  a2  a3  Ms  ~  P 


(4-74) 


b4  =  aj  ag  M< 


(4-75) 


Equation  4-80  can  be  decomposed  into  real  and  imaginary  parts. 


and  the  solution  for  c'  is  of  the  form 


where  dj  and  d3  are  related  to  d3  and  d4  as  follows! 


(4-76) 


(4-77) 


K  '  C92kS 


(4-78) 


84  =  1  CgKc 

The  solution  of  Equation  4-65  is  given  by 


(4-79) 


m2  =  b2  i  (b*  -  4a4)1/2 


(4-80) 


(c*)2  =  d3  +  id4 


(4-81) 


c'  =  dj  +  i d  2 


(4-82) 


dj  ■  ~-[(d32  +  d42)Jf  +  d3]* 


(4-83) 
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(4-84) 


d2  =  -  j^=r-[(d32  +  d42)^  -  d3]^ 

Solving  Equations  4-51  and  4-82  for  the  wavespeed,  c,  and  the  damping 
constant,  af  yields  the  following: 


(4-85) 


a 


(4-86) 


Because  of  the  nature  of  the  quadratic  formula  in  Equation  4-80,  two 
solutions  are  generated  for  each  set  of  parameters.  These  two  solutions  are 
called  waves  of  the  first  kind  and  waves  of  the  second  kind.  The  following 
subsection  presents  two  verifications  of  the  above  equation. 


4.3.3  Wave  Propagation  Velocity  for  Special  Cases 

In  this  subsection  plane  wave  propagation  velocities  are  determined 
for  three  special  cases: 

•  incompressible  grains  with  generalized  Darcy  flow  (r  «  0) 
and  no  damping  (a  ■  0); 

•  the  product  of  the  permeability  and  the  excitation  frequency 
approaches  zero; 

•  the  product  of  the  permeability  and  the  excitation  frequency 
approaches  infinity. 


In  the  latter  two  cases,  it  will  be  demonstrated  that  the  effect  of  per- 
meablity  and  excititation  frequency  on  the  wave  velocity  are  equivalent  and 
that  low  values  of  the  product  give  the  lower  bound  value  of  the  wavespeed  and 
that  high  values  give  the  upper  bound  wavespeed. 
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4, 3. 3.1  Incompressible  Grains,  Generalized  Darcy  Flow 


For  infinitely  stiff  grains.  Kg  -*  »  or  Cg  ®  0  and  for  generalized 
Oarcy  flow  the  mass  increment  factor,  r,  equals  zero.  The  fluid  friction  is 
simply  proportional  to  the  relative  fluid  velocity  and  is  independent  of  the 
excitation  frequency.  Equations  4-76  through  4-79  can  be  rewritten  as 


« 


n 


83 


s 


1_ 


«4  »  1 

Equations  4-68  through  4-71  yield 

pf  pg 

bi  «  (l  -  n) 

np 


(1-n)2  K,  Pf  *  n(l-n)  Kf  Pg  +  n  Pf  N$ 

b3 - ITkJ - 


Kf  ♦  nMs 


(4-87) 

(4-88) 

(4-89) 

(4-90) 

(4-91) 

(4-92) 

(4-93) 

(4-94) 


66 


Substitution  of  Equations  4-91  and  4-92  into  Equation  4-66  gives 


(4-95; 

Substitution  of  Equations  4-93  and  4-94  into  Equation  4-67  gives 

(4-96) 

Equations  4-95  and  4-96  are  identical  to  equations  derived  by  van  der  Kogel 
(1977)  and  serves  as  a  check  on  an  elementary  application  of  our  equations. 
Substitution  of  these  equations  into  Equation  4-80  and  4-85  with  a  =  0  gives  the 
wavespeeds  of  the  first  and  second  kinds  for  incompressible  grains  and  genera¬ 
lized  Oarcy  flow. 

4. 3. 3. 2  Lower  Bound  Wavespeed 

Whenever  the  permeability  and/or  the  excitation  frequency  approach 
zero,  there  is  no  relative  motion  between  the  pore  fluid  and  the  porous  skele¬ 
ton  and  an  undrained  loading  situation  exists.  Substitution  of  ku  ■  0  into 
Equations  4-76  through  4-79  yields 

aj  ■  oo  (4-97) 

Pf 

a2  •  ~  (4-98) 


+  (i-h).  _  !$s_ 

"o  Kg2 


84  ■  1 


(4-99) 

(4-100) 


Substitution  of  Equations  4-97  through  4-100  into  Equation  4-68  through  4-71 
result  in  no  change  in  the  expressions  for  bj  and  83,  but  b2  and  b4  are  given 
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by 

b3  =  »  (4-101) 

b4  =  oo  (4-102) 

Substitution  of  bi  through  b4  into  Equations  4-66,  '-67,  4-65  and  4-85  gives 


c2 


,  a42  +  a3  Hs 

cz  -  - 

P  *3 


0 


(4-103) 


Solution  of  Equation  4-103  yields 


c2  ■  0 


(4-104) 


and 


c2 . 

P®3 


(4-106) 


Equation  4-104  indicates  that  propagation  velocity  of  waves  of  the  second  kind 
equals  zero  for  the  undrained  case.  Substitution  of  Equations  4-99  and  4-100 
into  Equation  4-105  yields 


c2‘  (p)Hf 


(4-106) 


where  Mf  is  the  fully  coupled  undrained  constrained  modulus  which  is  given  by 


Hf  ■*•♦«$- 


Kg 


♦  K»  *5 


Ka  Ks 

K« - xr~  -  Kf 


K*  K® 


(4-107) 


68 


This  is  identical  to  the  fully  coupled  constrained  modulus  derived  by  Blouin 
and  Kim  (1984)  for  undrained  uniaxial  compression  (see  Section  2,  Equation 
2-17). 


4. 3. 3. 3  Upper  Bound  Wavespeed 

Whenever  the  permeability  and/or  excitation  frequency  approach  infi¬ 
nity,  Equations  4-69,  4-71  and  4-76  are  given  by 

b2  =  0  (4-108) 


b4  =  0  (4-109) 

aj  =  0  (4-110) 

The  Equations  for  bj,  03,  a2,  83  and  a4  remain  unchanged.  Substitution  of 
these  expressions  for  a's  and  b's  into  Equation  4-66  and  4-85  gives 


(4-111) 


where  the  plus  sign  gives  the  velocity  of  waves  of  the  first  kind  and  the 
minus  sign  gives  the  velocity  of  waves  of  the  second  kind.  The  constants  bj 
and  b3  in  Equation  4-111  are  given  by 


and 


*1 


n_  lH|j 
*f  *g  "  Kg2 


(4-112) 


P  -  2pf 


"•te 


b3 


(4-113) 


4.4  Implementation  Into  Computational  Algorithm 


In  order  to  fully  utilize  the  analytic  expression  for  wave  velocity  and 
damping  derived  in  the  previous  subsection,  the  general  solutions  for  the 
wavespeeds  and  damping  given  by  Equations  4-85  and  4-86  have  been  coded  into  th< 
program  TWAVE,  a  listing  of  which  is  included  in  Appendix  C.  In  Section  5 
TWAVE  is  used  to  evaluate  the  influence  of  various  parametric  changes  in 
material  properties  and  excitation  frequencies  on  wavespeed  in  two  phase 
media. 

4.4.1  Computational  Algorithm  (TWAVE) 

The  TWAVE  imput  is  contained  in  an  input  file  named  "TAPES" 
consisting  of  five  "cards."  As  shown  in  Tables  4.1  and  4.2  the  cards  are 
arranged  as  follows: 

Card  1,  frequency  option; 

Card  2,  excitation  frequency; 

Card  2,  pore  properties  including  porosity,  permeability, 
mass  increment  factor,  flow  path  constant; 

Card  4,  mass  densities  of  the  pore  fluids  and  solid  grains; 

Card  5,  bulk  moduli  of  the  pore  fluid,  solid  grains  and  porous 
skeleton;  constrained  modulus  of  the  porous  skeleton. 

The  input  values  are  in  floating  point  format  within  the  columns  shown  in 
Table  4.1.  Any  consistent  set  of  units  can  be  utilized. 

Host  of  the  input  parameters  listed  in  Table  4.2  are  straightforward 
except  for  the  frequency  option,  NF,  the  mass  increment  factor,  r,  and  the 
flow  path  constant,  r. 

•  NF  For  any  given  set  of  properties  the  speed  and  damping  of 
waves  of  the  first  kind  vary  as  a  function  of  frequency. 
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At  the  frequency  where  the  wavespeed  has  the  greatest 
variation,  damping  is  at  a  maximum.  The  program  computes 
and  uses  this  approximate  critical  frequency  and  computes 
the  corresponding  wavespeed  and  damping  when  NK  =  1. 

♦  r  The  mass  increment  factor,  as  explained  by  Kim  and  Blouin,  1984, 
is  a  frequency  dependent  flow  resistance  term  resulting  from 
fluid  friction.  It  is  treated  as  an  additional  inertial 
mass  term.  For  the  simple  extension  of  Darcy  flow,  r  equals 
zero,  which  is  termed  generalized  Darcy  flow.  In  Blot's  more 
rigorous  theoretical  treatment,  based  on  laminar  flow  conditions, 
mass  increment  factors  are  introduced  which  are  dependent  on 
the  pore  geometry.  For  idealized  pore  geometries  characterized 
by  circular  and  flat  ducts,  the  mass  increment  factors  have  values 
of  1/3  and  1/5,  respectively.  Whenever  the  approximation  using 
r  is  utilized,  the  flow  path  constant,  r,  must  be  specified  as 
zero. 

The  flow  path  constant,  r,  represents  Biot's  exact  expression 
for  frequency  dependent  fluid  friction  in  laminar  flow.  As 
with  the  approxiamate  solution,  values  for  Darcy  flow  and 
circular  and  flat  ducts  are  specified.  These  values  are  0,  1 
and  firi,  respectively.  Whenever  the  exact  solution  using 
r  is  utilized,  the  approximation,  r,  must  be  set  to  zero. 
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Table  4.1.  Summary  of  input  parameters  and  FORTRAN  equivalent. 


Summary  of  Input  Parameters 


FORTRAN  Equivalent 


Card  1 


Card  2 


Card  3 


Card  4 


Table  4.2.  Description  of  input  parameters. 


-  Frequency  Options 
NF  =  0  Use  given  input  frequency 

NF  =  1  Use  internally  computed  approximate  critical  frequency 
(neglect  input  frequency) 


-  Loading 

w  =  Frequency  { rad/sec) 

-  Pore  Properties 

n  =  Porosity 

k  =  Permeability  coefficient  (ex.  in/sec) 

r  *  Mass  increment  factor  (from  Biot's  1961  paper) 

r  *  0  (generalized  Darcy  flow) 

r  •  1/5  (flat  pore) 

r  »  1/3  (circular  pore) 

r  *  Flow  path  constant 

r  *  0 _  (generalized  Darcy  flow) 

r  •  /m  (flat  pore) 

r  ■  1  (circular  pore) 

-  Hass  Densities 

yw  •  Unit  weight  of  standard  water  (ex.  62.4/(123)  lb/in3) 
g  *  Gravitational  acceleration  (ex.  (32.2) (12)  in/sec?) 

.G. f  o  Specific  gravity  pore  fluid 
i.G.g  *  Specific  gravity  of  solid  grain 
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Table  4.2.  (concluded) 


Card  5  -  Moduli 

Kf  =  Bulk  modulus  of  pore  fluid 

Kg  =  Bulk  modulus  of  solid  grain 

Ks  =  Bulk  modulus  of  porous  skeleton 

Ms  =  Constrained  modulus  of  porous  skeleton 
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Conservation  of  Fluid  Hass 
n  p.  Vt  a  n’  p*'  vt' 

Vt  =  apparent  fluid  volume  before  compression 

Vt'  =  (1  +  ep)  Vt:  apparent  fluid  volume  after  compression 

ev  =  volumetric  strain  of  porous  skeleton 

€p  =  volumetric  diffusion  of  pore  fluid 


Vt 


I 


Before  Compression 


After  Compression 


figure  4,2.  Schematic  illustration  of  conservation  of  pore 
fluid  mass  in  saturated  porous  materials. 
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SECTION  5 


NUMERICAL  ANALYSIS  OF  WAVE  PROPAGATION 
IN  SATURATED  POROUS  MEDIA 


5.1  INTRODUCTION 

In  this  section  the  analytic  closed  form  solution  for  wave  velocity  and 
damping  in  fully  saturated  elastic  porous  media,  as  implemented  in  TWAVE,  is 
used  to  evaluate  the  influence  of  excitation  frequency  and  systematic 
variations  in  material  properties  on  the  propagation  velocity  and  damping. 
Important  aspects  of  this  parameter  study  are  investigations  of 

•  velocities  and  damping  for  waves  of  both  the  first  and  second 
kind  and  variation  of  velocity  and  damping  as  functions  of  the 
permeabi 1 i ty-f requency  product ; 

•  the  influence  of  pore  shape  and  fluid  friction  on  the  wave  prop¬ 
agation  characteristics; 

•  The  influence  of  water  properties  (sea  water  vs  fresh  water), 
porosity,  and  skeleton  properties  (Poisson's  ratio,  bulk  modulus) 
on  the  wave  propagation  characteristics;  and 

•  the  adequacy  of  the  approximate  expression  for  fluid  friction 
(used  in  the  numerical  model)  as  compared  to  the  exact  expression 
for  fluid  friction. 

The  overall  results  of  this  study  provide  the  basis  for  quantitative  eva¬ 
luation*  and  judgments  regarding  wave  velocity  characteristics  and  spatial 
attenuation  in  a  broad  variety  of  real  earth  materials.  They  are  valuable 
guides  to  our  basic  understanding  of  wave  propagation  phenomena  and  to 
planning  and  analysis  of  the  future  numerical  studies  in  this  program. 
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5.2  INFLUENCE  OF  POROSITY  AND  SKELETON  MODULUS  ON  WAVESPEED  AND  DAMPING 

As  described  in  Section  4,  the  expression  of  Equation  4-82  for  the 
wavespeed  in  saturated  porous  media  generates  two  solutions  which  are  referred 
to  as  waves  of  the  first  and  second  kind.  The  uave  of  the  first  kind  is  the 
common  compressional  wave  through  the  porous  skeleton-fluid  matrix  with 
stresses  distributed  within  the  matrix  in  proportion  to  the  properties  of  the 
components.  The  wave  of  the  second  kind  is  not  well  understood,  but  appears 
to  be  associated  with  a  surge  of  pore  fluid  through  the  matrix.  We  plan  on 
investigating  the  phenomena  associated  with  waves  of  the  second  kind  with  the 
help  of  numerical  simulations  during  the  coming  year. 

In  this  subsection,  TWAVE  is  used  to  generate  plots  of  wavespeed  and 
damping  as  a  function  of  the  frequency-permeability  product  for  a  variety  of 
material  properties  chosen  to  envelop  the  properties  of  most  porous  saturated 
rocks  and  soils.  The  range  of  these  parametric  values  are  listed  in  Table  5.1. 

The  results  of  the  parametric  calculations  are  presented  in  the  form  of 
wavespeed  and  damping  as  a  function  of  the  product  of  frequency  and  per¬ 
meability.  Analysis  of  the  solutions  presented  in  Section  4  shows  that  the 
wavespeed  and  specific  damping  are  a  function  of  the  frequency  permeability 
product.  As  long  as  the  product  of  the  two  doesn't  change,  there  will  be  no 
change  in  either  specific  damping  or  wavespeed.  Thus  the  product  of  the  two 
parameters  makes  a  convenient  plotting  reference. 

The  damping  coefficient,  a,  from  Equations  4-47  through  4-50  is  nor¬ 
malized  by  the  excitation  frequency,  w,  and  is  termed  the  specific  damping. 

It  is  a  measure  of  the  rate  of  motion  and  stress  attenuation  per  wavelength 
with  increasing  depth. 

A  set  of  parametric  calculations  using  TWAVE  was  run  on  porous  material 
saturated  with  fresh  water.  Skeleton  modulus  was  varied  by  order  of  magnitude 
increments  form  1,000  psi  to  1,000,000  psi.  Poisson's  ratio  of  the  skeleton 
was  fixed  at  0.3.  Porosities  of  0.2,  0.4  and  0.6  were  computed.  Values  of 
flow  path  constant,  r,  for  the  exact  solution  of  0,  fljz,  1  and  2  were  used. 
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The  numerical  data  are  organized  in  three  groups,  as  shown  in  Table  5.2, 
according  to  the  three  different  porosities.  For  each  porosity  four  sets  of 
wavespeed  and  damping  plots  are  presented,  one  for  each  of  the  four  skeleton 
moduli.  Each  Figure  has  plots  of  the  four  flow  paths. 

5.2.1  Speed  of  Waves  of  the  First  Kind 

•  There  is  a  wavespeed  transition  from  lower  to  higher  speed  with 
increasing  wk,  The  lower  bound  wavespeed,  as  discussed  in 
Section  4. 3. 3. 2,  is  identical  to  the  wavespeed  computed  from 
the  undrained  constrained  modulus  using  Equation  3-12.  At  the 
lower  bound  values  of  wk  there  is  no  significant  relative  motion 
between  the  pore  water  and  the  porous  skeleton.  Wavespeeds 

at  higher  values  of  wk  approach  an  upper  bound  discussed  in 
Section  4. 3. 3. 3.  The  transition  occurs  over  the  50  to  10,000 
range,  which  for  a  typical  sand  porosity  of  0.1  in/s  corresponds 
to  a  frequency  range  of  about  8  to  1,600  cycles  per  second. 

•  The  wavespeed  transition  occurs  more  sharply  for  the  pores  having 
less  fluid  friction  (minimum  r). 

•  The  wavespeed  increases  with  increasing  skeleton  modulus  but  the 
magnitudes  of  the  difference  between  the  upper  and  lower  bound 
wavespeeds  do  not  vary  consistently  with  increasing  modulus. 

•  There  appears  to  be  little  influence  of  porosity  on  the  wavespeed 
increase  for  soft  skeletons;  but  there  is  a  large  jump  in  wavespeed 
for  high  porosity  materials  with  stiff  skeletons  and  almost  no 
jump  in  wavespeed  for  low  porosity  stiff  materials. 

5.2.2  Speed  of  Waves  of  the  Second  Kind 

•  There  is  a  wavespeed  transition  from  near  zero  to  an  upper  bound 
velocity  with  increasing  values  of  wk.  The  lower  bound  approaches 
assymptotically  to  zero. 
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•  The  transition  occurs  more  sharply  for  the  pores  having  lower  friction 
resistance,  and  it  occurs  at  a  somewhat  lower  uk  than  the  transition 
for  the  waves  of  the  first  kind. 

•  The  magnitude  of  the  wavespeed  jump  increases  with  increasing  skeleton 
stiffness. 

•  There  is  a  relatively  small  influence  of  porosity  on  the  wavespeed 
increase,  compared  to  a  strong  influence  of  skeleton  stiffness. 

5.2.3  Damping  of  Waves  of  the  First  Kind 

•  Specific  damping  shows  a  strong  peak  at  the  point  of  steepest  increase 
in  wavespeed. 

•  There  is  no  damping  at  low  values  of  uk,  where  there  is  no  restive 
motion  between  the  pore  water  and  skeleton.  At  higher  values  of  uk 
the  damping  appears  to  be  decreasing  toward  zero. 

•  At  the  damping  peak  the  pores  with  less  frictional  resistance  have 
higher  specific  damping  than  those  with  higher  frictional  resis¬ 
tance.  At  higher  values  of  uk  this  trend  is  reversed,  and  the 

higher  frictional  pores  have  the  most  damping.  At  low  values  of  uk  the 
damping  in  all  cases  is  about  the  same. 

•  The  specific  damping  is  strongly  influenced  by  the  skeleton  modulus, 
but  as  explained  in  Section  5.6.2  the  damping  drops  to  zero  at 
certain  modulus-porosity  combinations. 

•  Porosity  appears  to  have  little  influence  on  specific  damping. 

5.2.4  Damping  of  Waves  of  the  Second  Kind 

•  Specific  damping  of  waves  of  the  second  kind  is  extremely  strong  at  low 
values  of  uk  and  drops  exponentially  toward  zero  at  higher  uk  values. 
Damping  is  near  zero  at  higher  uk  values.  Damping  is  near  zero  once  the 
wavespeeds  approach  their  limiting  upper  bound,  but  as  the  wavespeeds 


drop  below  the  transition  zone  the  damping  becomes  very  strong.  Thus, 
we  speculate  that  it  may  be  very  difficult  to  detect  waves  of  the 
second  kind  below  their  upper  bound  velocities.  Note  that  the  magni¬ 
tudes  of  the  specific  damping  plots  are  about  three  orders  of  magnitude 
higher  than  the  corresponding  plots  for  damping  of  waves  of  the  first 
kind. 

•  There  is  almost  no  influence  of  the  flow  path  factor,  r,  on  the 
specific  damping. 

•  Increasing  skeleton  modulus  significantly  decreases  the  specific 
damping. 

•  Porosity  has  very  little  influence  on  the  specific  damping. 

5.3  COMPARISON  OF  FRESH  WATER  AND  SEA  WATER 


Since  many  porous  materials  are  saturated  with  sea  water,  wave  velocities 
and  damping  for  a  typical  sand  saturated  with  sea  water  were  computed  for  com¬ 
parison  to  similar  material  saturated  with  fresh  water.  For  sand  having  pro¬ 
perties  listed  in  Table  5.2,  with  n  »  0.4  and  ks  »  10,000  psi,  the  series  of 
computations  shown  in  Figure  5.13  was  generated  using  a  bulk  modulus  for  sea 
water  of  350,000  psi.  Comparison  of  these  results  to  the  corresponding  plots 
for  fresh  water  in  Figure  5.6  shows  that: 


•  wavespeeds  for  waves  of  the  first  kind  are  higher  in  the  sea  water 
saturated  material  because  of  the  higher  bulk  modulus  of  the  sea 
water  and  the  correspondingly  higher  bulk  modulus  of  the  mixture? 

•  wavespeeds  for  waves  of  the  second  kind  are  not  noticeably  affected 
by  the  modulus  of  the  pore  water? 

•  specific  damping  of  the  first  kind  is  lower  in  the  sea  water 
saturated  material? 

•  however,  specific  damping  for  waves  of  the  second  kind  is  higher 
in  the  sea  water  saturated  material. 
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The  reasons  for  the  last  three  observations  are  not  well  understood  at  this 
point. 

5.4  INFLUENCE  ON  POISSON'S  RATIO  ON  VELOCITY  AND  DAMPING 

The  influence  of  Poisson's  ratio  of  the  skeleton  was  determined  by  com¬ 
puting  the  wavespeeds  and  damping  for  a  typical  saturated  sand  having  the  pro¬ 
perties  listed  in  Table  5.2  for  n  =  0.4  and  ks  =  10,000  psi,  but  with 
Poisson's  ratios  of  0.2  and  0.4.  As  Equation  2-25  demonstrates,  raising 
Poisson's  ratio  is  equivalent  to  reducing  the  constrained  skeleton  modulus. 

The  results,  plotted  in  Figure  5.14  for  v  =  0.2  and  5.15  for  v  =  0.4,  can  be 
compared  to  Figure  5.6  for  p  =  0.3.  Variation  in  Poisson's  ratio  has  only  a 
slight  effect  on  both  velocity  and  specific  damping. 

•  Wavespeed  for  both  waves  of  the  first  and  second  kind  are  slightly 
higher  with  lower  Poisson's  ratio. 

•  Specific  damping  for  both  wave  types  is  slightly  lower  with  lower 
Poisson's  ratio. 

5.5  COMPARISONS  BETWEEN  EXACT  AND  APPROXIMATE  EXPRESSIONS  FOR  FLUID  FRICTION 

As  explained  in  Section  4.4.1,  the  numerical  codes  TPDAP  II  and  MPDAP  use 
an  approximation  of  r  the  frequency  dependent  fluid  friction  given  by  the  mass 
increment  factor,  (Kim  and  Blouin,  1984).  This  approximation  is  necessary 
because  the  fluid  friction  term  is  frequency  dependent  and  the  frequency  is 
not  known  beforehand.  Hence,  either  an  iterative  scheme  or  an  approximation 
such  as  the  above  must  be  employed.  An  iterative  scheme  would  be  too  time 
consuming  for  most  applications. 

In  order  to  check  the  adequacy  of  the  mass  increment  approximation  for 
fluid  friction,  we  performed  a  suite  of  calculations  on  "standard"  saturated 
material  with  properties  shown  in  Table  5.2  with  a  porosity,  n,  of  0.4.  The 
results  are  plotted  in  Figure  5.16  through  5.19  for  values  of  skeleton  modulus 
ranging  from  1,000  psi  to  1,000,000  psi.  Each  plot  compares  the  exact 
expression  for  fluid  friction  in  a  circular  duct  (r  ®  1)  with  the  approximate 
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exp  e-osion  (r  =  1/3).  These  are  also  compared  to  the  generalized  Oarcy  fluid 
fric- ion  approximation  (r  =  r  =  0).  Conclusions  determined  from  analysis  of 
these  plots  include: 

®  The  mass  increment  approximation  gives  an  excellent  representation  of 
wavespeed  and  damping  up  to  a  criteral  value  of  wk,  where  the  damping 
and  the  change  in  wavespeed  are  a  maximum.  At  higher  values  of  the 
frequency  permeability  product,  the  mass  increment  approximation  becomes 
increasingly  less  accurate. 

•  For  wk  products  above  the  critical  value,  the  mass  increment  approxima¬ 
tion  underestimates  the  wavespeeds  of  both  kinds  of  waves  and  under¬ 
estimates  the  damping  for  waves  of  the  first  kind.  Damping  of  waves 
of  the  second  kind  appears  to  be  independent  of  the  fluid  friction 
term.  The  upper  bound  wavespeed  is  independent  of  the  flow  path  con¬ 
stant,  r,  and  the  wavespeeds  converge  to  the  Darcy  flow  case. 

However,  in  the  approximation,  the  mass  increment  factor,  as  shown  in 
Equation  4-19,  is  still  carried  at  higher  values  of  wk. 

5.6.  INFLUENCE  OF  POROSITY  AND  SKELETON  MODULUS  ON  DAMPING  AND  WAVESPEED 

Section  5.2  described  the  wavespeeds  and  specific  damping  as  a  function 
of  the  frequency-permeability  product.  In  this  section,  we  evaluate  the 
wavespeeds  and  damping  as  a  function  of  the  skeleton  bulk  modulus.  Material 
properties  used  for  this  study  are  the  same  as  those  listed  in  Table  5.2.  For 
each  porosity  we  evaluate  the  upper  and  lower  bound  wavespeeds  which  are  inde¬ 
pendent  of  the  flew  path  constants.  In  addition,  we  evaluate  the  values  of 
the  critical  frequency-permeability  product  where  the  specific  damping  of 
waves  of  the  first  kind  reaches  a  the  maximum.  At  this  critical  frequency- 
permeability,  we  compute  velocities  of  waves  of  the  first  kind  and  specific 
damping  of  waves  of  both  the  first  and  second  kind.  Figures  5.20  through  5.22 
show  the  results  for  n  =  0.2,  n  »  0.4  and  n  ■  0.6,  respectively.  Each  figure 
includes  plots  for  the  four  flow  path  constants  listed  in  Table  5.2. 
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5.6.1  Critical  Frequency-Permeabi 1 ity 

•  Critical  wk  shows  a  peak  at  higher  skeleton  bulk  moduli 

(Ks  >  100,000  psi).  The  value  -of  Ks  where  the  peak  occurs  decreases 
as  the  porosity  increases. 

•  The  magnitudes  of  critical  wk  increase  with  increasing  porosity. 

•  However,  the  critical  wk  values  do  not  vary  monitonically  with 
variations  in  fluid  friction.  The  lowest  value  of  wk  is  for  an 
r  value  of  1.0;  values  both  above  and  below  1.0  produce  higher 
values  of  the  critial  wk. 

5.6.2  Maximum  Specific  Damping  of  Waves  of  the  First  Kind 

•  The  specific  damping  curve  has  a  saddle  point,  reaching  zero 
damping  at  the  skeleton  bulk  modulus  where  the  critical  wk 
reaches  its  peak. 

•  The  magnitude  of  specific  damping  increases  with  increasing 
porosity. 

•  The  magnitude  of  specific  damping  increases  with  decreasing 
fluid  friction  (decreasing  r). 

5.6.3  Specific  Damping  of  Waves  of  the  Second  Kind  at  Critical  wk. 

•  Specific  damping  at  critical  wk  decreases  exponentially  with 
increasing  skeleton  bulk  modulus. 

•  The  magnitude  of  specific  damping  decreases  as  the  porosity  increases. 

•  The  magnitude  of  specific  damping  is  lowest  for  material  where 
the  flow  path  constant,  r,  is  equal  to  2.0. 

5.6.4  Velocities  of  Waves  of  the  First  Kind 

•  Up  to  a  skeleton  modulus  Ks,  of  about  100,000  psi,  wave  velocities 
increase  gradually  with  increasing  modulus.  Velocities  are 
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higher  in  the  materials  having  lower  porosity.  Around  the  value 
of  Ks  where  the  critical  «k  reaches  its  peak,  wave  velocities 
are  increasing  exponentially  and  converge  to  a  single  value  where 
critical  wk  is  a  maximum. 

•  Both  maximum  and  minimum  wave  velocities  are  the  same  at  the  point 
where  the  damping  of  waves  of  the  first  kind  drops  to  zero. 

•  For  the  flow  path  constant  r  =  2,  the  wave  velocities  at  critical 
ok  are  approximately  the  average  of  the  maximum  and  minimum  values. 

For  the  other  values  of  flow  path  constants  (r  »  0,  /2/3,  and  1), 

the  wave  velocities  at  critical  tok  are  slightly  lower  than  the 
average  of  the  maximum  and  minimum  values, 

5.6.5  Maximum  Velocities  of  Waves  of  the  Second  Kind 

•  Wave  velocities,  in  general,  increase  as  the  skeleton  bulk  modulus 
increases  and  the  magnitudes  of  the  wave  velocities  are  higher  in  the 
materials  having  higher  porosity. 

•  The  maximum  wave  velocity  is  not  influenced  by  the  fluid  friction. 

5.7  TPOAP  II  NUMERICAL  CALCULATIONS 

A  series  of  3  one-dimensional  calculations  of  a  vertically  propagating 
planar  compression  wave  were  performed  using  the  two  phase  finite  element 
program  TPDAP  II.  The  input  loading,  as  shown  in  Figure  5.23,  was  a  short 
rise  time  triangular  pulse  with  a  peak  stress  of  5,000  psi  and  a  positive 
phase  duration  of  10  msec.  The  loading  pulse  was  applied  to  saturated  sand 
having  the  properties  listed  in  Figure  5.23.  Snapshots  of  the  pore  water 
pressure  profiles  at  four  different  times,  from  10  to  40  msec,  are  shown  in 
Figure  5.24.  Three  calculations  are  shown  for  assumed  permeabilities  of 
0.001,  0.1  and  10  in/s. 

The  trends  in  these  calculations  are  similar  to  those  observed  in  the 
parametric  calculations  of  the  previous  subsections.  The  wavespeed  increases 
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substantially  with  increasing  permeability.  For  the  lowest  permeability  the 
velocity  of  the  wavefront  is  about  5300  ft/s  and  for  the  highest  permeability 
it  is  about  5900  ft/s.  The  wavespeed  computed  for  this  material  according  to 
the  decoupled  undrained  modulus  described  by  Blouin  and  Kim  (1984)  is  about 
5200  ft/s  which  agrees  very  well  with  the  low  permeability  calculation. 

There  is  a  dramatic  alteration  of  the  wave  shape  for  the  intermediate 
permeability  calculation.  The  wavefront  is  smeared  and  the  amplitude  is 
significantly  attenuated  relative  to  the  calculations  in  the  less  permeable 
and  more  permeable  materials.  As  discussed  in  the  previous  subsections, 
excess  damping  occurs  in  the  transition  region  between  the  lower  bound  and 
upper  bound  wavespeeds. 
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Table  5.1.  Material  properties  used  in  parameter  study. 


Fixed  Constants: 

Grain 

Specific 

Gravity 

Bulk 

Gravity  of 

Density 

Constant 

Modulus 

Grains 

yw  =  62.4  lb/ft3, 

g  =  32.17  ft/sec2. 

Kg  =  5  x  10®  psi. 

S.G.g  =  2.6’ 

Flow  Path  Constant  (exact  sxpression  for  fluid  friction): 
r  =  0,  fijz,  1.0,  2.0 

Mass  Increment  Factor  (approximate  expression  for  fluid  friction): 
r  ■  0,  1/5,  1/3,  and  1 

Properties  of  Pore  Fluid: 


Fluid 

Specific  Gravity 

Bulk  Modulus 

of  Pore  Fluid 

Fresh  Water  - 

Kf  •  0.29  x  10®  psi 

and 

S.G.f  ■  1.0 

Sea  Wator 

Kf  *  0.35  X  10®  psi 

and 

S.G.f  *  1.026 

Properties  and  Porous  Skeleton: 

Bulk  Modulus  -  Kg  **  1,000,  10,000,  100,000,  1,000,000  psi 

Poisson's  ratio  -  v  =  0,2,  0.3,  0.4 

Porosities: 

n  ■  0.2,  0.4,  and  0.6 
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Table  5.2.  List  of  figures  for  influence  of  porosity  and  skeleton  modulus 
on  wavespeed  and  damping. 


Skeleton  Modulus  : 


n  =  0.2 

n  a  0.4 

n  »  0.5 

Kg  =  1,000 

Fig.  5.1 

Fig.  5.5 

Fig.  5.9 

Kg  =  10,000 

Fig.  5.2 

Fig.  5.6 

Fig.  5.10 

Kg  =  100,000 

Fig.  5.3 

Fig.  5.7 

Fig.  5.11 

Kg  *  1,000,000 

Fig.  5.4 

Fig.  5.8 

Fig,  5.12 

Constants ; 

Kf  *  0.29  x  106  psi  (Bulk  Modulus  of  Pore  Water) 
v  a  0.3  (Poisson's  ratio) 
r  a  o,  ✓573 r  1,  2  (Flow  Path  Constant) 
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Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Fneq-Penm  Product 


Figure  5.1a.  Wave  velocities,  n  =  0.2  and  K  °  1,000  psi 
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Specific  Damping  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 
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Figure  5.1b.  Specific  damping,  n  *  0.2  and  K$  =  1,000  psi. 
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Figure  5.3a.  Wave  velocities,  n  *  0.2  and  Ks  ■  100,000  psi. 
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-ALPHA2/W  [(1/in) /(rad/s)]  “ALPHA1/w  [  (1/in)  /  (rad/s)  ] 


Specific  Damping  of  First  and  Second  Kind 
as  a  Function  of  Freq-Penm  Product 


0.1  1  10  100  1000  10000 


w#k  [  (rad/s)  #  (in/s)  ] 

Figure  5.3b.  Specific  damping,  n  6  0.2  and  K  =  100,000  psi. 
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Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


Figure  5.4a.  Wave  velocities,  n  9  0.2  and  K  *  1,000,000  psi 
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Specific  Damping  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 
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Figure  5.4b.  Specific  damping,  n  a  0.2  and  Ks  °  1,000,000  psi. 
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Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


Figure  5.5a.  Wave  velocities,  n  =  0.4  and  K  *  1,000  psi 


Specific  Damping  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


Figure  5.5b.  Specific  damping,  n  a  0.4  and  K  =  1,000  psi 
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Figure  5.6a.  Wave  velocities,  n  *  0.4  and  K  °  10,000  psl 
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Specific  Damping  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


w*k  [  (rad/s)  *  (in/s)  ] 


Figure  5.6b.  Specific  damping,  n  =  0.4  and  K$  *  10,000  psi. 
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Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


Figure  5.7a.  Wave  velocities,  n  *  0.4  and  K  =  100,000  psi 
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Figure  5.8b.  Specific  damping,  n  =  0.4  and  K  =  1,000,000  psi. 
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[(ft/s)] 


Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 
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Figure  5.9a.  Wave  velocities  In  porous  material  saturated  with  fresh 
water;  n  ■  0.6,  K,  °  1,000  psl. 
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Figure  5.9b.  Specific  damping  *n  porous  material  saturated  with 
fresh  water,  n  «  0.6,  K  •=  1,000  psi. 
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Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Penm  Product 


5200 


^  5000 


w*k  [  (rad/s)  #  (in/s)  ] 

Figure  5.10a.  Wave  velocities  in  porous  material  saturated  with 
fresh  water,  n  *  0.5,  K,  »  10,000  psl. 
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Figure  5.10b.  Specific  damping  in  porous  material  saturated  with 
fresh  water,  n  =  0.6,  K  -  10,000  psi. 
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C2  [(ft/s)]  Cl  [(ft/s) 


Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


w#k  [  (rad/s)  #  (in/s)  ] 

Figure  5.11a.  Wave  velocities  in  porous  material  saturated  with  fresh 
water,  n  *  0.6,  K$  *  100,000  psi. 
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-ALPHA2/w  [  (i/in)  /  (rad/s)  ] 


Specific  Damping  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 
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Figure  5.11b,  Specific  damping  in  porous  material  saturated  with 
fresh  water,  n  *  0.6,  K  =  100,000  psi. 
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Cl  t  (ft/s)  ] 


Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


w#k  [  (rad/s)  #  (in/s) ) 

\  . 

Figure  5.12a.  Wave  velocities  in  porous  material  saturated  with  fresh 
water,  n  ■  0.6,  Ks  *  1,000,000  psi. 
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Specific  Damping  of  First  and  Second  Kind 
as  a  Function  of  Fneq-Perm  Product 
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Figure  5.12b.  Specific  damping  in  porous  material  saturated  with 
fresh  water,  n  =  0.6,  K$  °  1,000,000  psi. 
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Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 
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Figure  5.13a.  Wave  velocities  In  porous  material  saturated  with 
sea  water,  n  *  0.4,  *  10,000  psi. 
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Figure  5.13b.  Specific  damping  in  porous  material  saturated 
with  sea  water,  n  »  0.4,  *  10,000  psi. 
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Specific  Damping  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


w#k  [  (rad/s)  #  (in/s)  ] 

Figure  5.14b.  Specific  damping  with  v  *  0.2  and  n  *  0.4, 
K  =  10,000  psi. 
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Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 
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figure  5.15b.  Specific  damping  with  v  *  0.4  and  n  «  0.4, 
K.  *  10,000  psi . 
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C2  [(ft/s))  C1  t  (ft/s)  ] 


Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


w*k  [  (rad/s)  #  (in/s)  ] 


Figure  5.16a.  Comparison  of  wavespeeds  computed  from  approximate  and 
exact  expressions  for  fluid  friction,  n  «  0.4, 

K  ■  1,000  psi. 
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Specific  Damping  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 
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Figure  5.16b.  Comparison  of  specific  damping  computed  from 
approximate  and  exact  expressions  for  fluid 
friction,  n  =  0.4,  K  =>  1,000  psi. 
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Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 


w#k  [  (rad/ s)  #  (in/s)  ] 

Figure  5.17a.  Comparison  of  wavespeeds  computed  from  approximate 
and  exact  expressions  for  fluid  friction,  n  *  0.4, 
K  ■  10,000  psi. 
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Figure  5.17b.  Comparison  of  specific  damping  computed  from 
approximate  and  exact  expressions  for  fluid 
friction,  n  B  0.4,  K  =  10,000  psi. 
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Wave  Veloc.  of  First  and  Second  Kind 


as  a  Function  of  Freq-Perm  Product 


w#k  [  (rad/s)  #  (in/s)  ] 

Figure  5.18a.  Comparison  of  wavespeeds  computed  from  approximate 
and  exact  expressions  for  fluid  friction,  n  ■  0.4, 
Kc  »  100,000  psl. 
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Figure  5.18b.  Comparison  of  specific  damping  computed  from  approximate 
and  exact  expressions  for  fluid  friction,  n  ■  0.4 
K  =  100,000  psi. 
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Wave  Veloc.  of  First  and  Second  Kind 
as  a  Function  of  Freq-Perm  Product 
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Figure  5.19a.  Comparison  of  wavespeeds  computed  from  approximate 
and  exact  expressions  for  fluid  friction,  n  *  0.4, 
K  *  1,000,000  psi. 
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Figure  5.19b.  Comparison  of  specific  damping  computed  from  approximate 
and  exact  expressions  for  fluid  friction,  n  *  0.4, 

K  =  1,000,000  psi. 
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Figure  5.20a.  Critical  values  of  to  k,  n  =  0.2. 


Maximum  Specific  Damping  of  First  Wave  as 
Function  of  Skeleton  Bulk  Modulus 
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Figure  5.20b.  Maximum  sped fie* damping  of  waves  of  the  first  kind,  n  =  0.2. 
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Figure  5.20c.  Specific  dasnping  of  waves  of  the  second  kind  at  critical  oj  k 
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Figure  5„2Ge.  Maximum  velocities  of  waves  of  the  second  kind 
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Figure  5.21b.  Maximum  specific  damping  of  waves  of  the  first  kind,  n  =  0.4. 
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Figure  5.21c.  Specific  damping  of  waves  of  the  second  kind  at  critical  w  k,  n  =  0.4. 
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Figure  5.21d.  Velocities  of  waves  of  the  first  kind,  n  =  0.04. 
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Figure  5.21e.  Maximum  velocities  of  waves  of  the  second  kind. 
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Figure  5.22a.  Critical  values  of  to  k,  n  =  0.6. 
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Figure  5.22b.  Maximum  specific  damping  of  waves  of  the  first  kind,  n  =  0.6. 
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Figure  5,22c.  Specific  damping  of  waves  of  the  second  kind  at  critical  w  k,  n  =  0.6. 
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Figure  5.22d.  Velocities  of  waves  of  the  first  kind,  n  =  0.04. 
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SECTION  6 


LIQUEFACTION  AND  CONSOLIDATION  UNDER  UNIAXIAL  STRAIN  LOADINGS 

6.1  INTRODUCTION 

It  is  well  established  that  both  near  surface  and  contained  explosions  in 
or  on  saturated  porous  soils  cause  extensive  zones  of  liquefaction  surrounding 
the  detonation  (Blouin  and  Shinn,  1983).  The  liquefaction,  in  turn,  dramati¬ 
cally  alters  the  late  time  ground  motions  and  cratering  processes.  Following 
the  onset  of  liquefaction,  ground  motions  in  the  vicinity  of  the  burst  change 
suddenly  from  high  frequency  motions  to  vefy  low  frequency  oscillatory 
motions  having  long  durations.  Immediately  following  formation  of  a  crater  or 
cavity  by  the  dynamics  of  the  explosion,  liquefied  material  flows  and/or 
slumps  in  to  partially  fill  the  excavation.  Following  the  gravity  controlled 
slumping  and  filling  process,  the  zone  of  liquefied  material  gradually  con¬ 
solidates,  with  excess  pore  water  working  its  way  out  to  the  ground  surface. 

The  consolidation  process  can  take  months,  depending  on  the  size  of  the  explo¬ 
sion  and  the  properties  of  the  in  situ  soil.  The  combined 
liquefaction/gravity  flow/consolidation  process  can  drastically  alter  the  size 
and  shape  of  craters  in  saturated  materials  and  has  been  a  factor  in  misin¬ 
terpretation  of  much  cratering  data.  In  the  past,  attempts  have  been  made  to 
relate  final  crater  shapes  and  volumes  directly  to  the  explosion  formation 
processes  without  regard  to  the  late-time  flow  and  consolidation  processes. 

In  order  to  gain  a  quantitative  understanding  and  validation  of  these 
processes  a  series  of  laboratory  tests  was  performed  which  attempted  to  repli¬ 
cate  important  aspects  of  the  in  situ  loading  conditions.  These  tests,  termod 
shock  consolidation  tests,  utilized  a  uniaxial  strain  load/unload/consolidation 
cycle  which  was  felt  to  be  appropriate  for  modeling  most  important  aspects  of 
the  in  situ  behavior.  Materials  cested  included  saturated  beach  sand, 
"undisturbed"  silt-sand-gravel  core  samples,  and  weakly  to  strongly  cemented 
porous  limestone  cores,  all  from  Enewetak  Atoll.  These  carbonate  materials 
were  furnished  by  the  Defense  Nuclear  Agency  from  core  taken  in  the  vicinity 
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of  the  thermonuclear  detonations  KOA  and  OAK.  Results  of  the  shock  con¬ 
solidation  tests  are  presented  in  this  section  and  are  compared  to  two  phase 
calculations  of  the  response  of  both  a  saturated  rock  and  soil  sample  prepared 
using  the  code  NKOCP  describea  in  Section  3. 

6.2  TEST  DESCRIPTION  AND  RESULTS 

A  series  of  17  shock  consolidation  tests,  listed  in  Table  6.1,  was  run  on 
Enewetak  beach  sand,  and  silt-sand-gravel  limestone  and  core  specimens.  As 
reported  by  Blouin,  Kim  and  Timian  (1986),  a  series  of  initial  tests  was 
attempted  in  a  triaxial  vessel,  but  as  soon  as  the  soil  liquefied  shortly  into 
the  unload  cycle,  it  became  impossible  to  control  the  unloading  and  con¬ 
solidation  portions  of  the  test.  Once  liquefaction  occurred,  the  sample 
behaved  as  a  fluid  and  could  support  only  hydrostatic  pressure.  During 
hydrostatic  unloading  it  was  impossible  to  control  the  radial  deformation  of 
the  sample  and  maintain  the  uniaxial  strain  condition. 

For  this  reason  all  testing  was  done  in  a  rigid  high  pressure  oedometer 
in  which  a  nominal  1  in  long  by  4  in  diameter  cylindrical  sample  was  loaded  in 
uniaxial  strain.  The  sample  was  constrained  laterally  by  the  cylindrical 
steel  walls  of  the  oedometer  vessel  and  was  loaded  from  the  top  by  means  of  a 
rigid  steel  piston.  High  pressure  seals  between  the  top  and  the  vessel  walls 
prevented  the  escape  of  pore  water.  Following  the  initial  undrained  load- 
unload  cycle,  drainage  was  provided  through  a  1/16  in.  line  in  the  bottom  of 
the  vessel.  Instrumentation  consisted  of  a  pore  pressure  gage  attached  to  the 
drainage  line,  a  linear  displacement  transducer  to  monitor  piston  displace¬ 
ment,  and  an  external  load  cell  to  measure  the  total  load  applied  to  the 
piston  by  the  load  frame.  A  schematic  section  view  of  the  oedometer  device  is 
shown  in  Figure  6.1. 

Important  steps  in  the  testing  procedure  include  the  following: 

1.  Sample  Saturation:  In  order  to  realistically  simulate  the  in  situ  con¬ 
ditions,  the  samples  must  be  completely  saturated.  If  complete  saturation 
isn't  achieved,  the  undrained  sample  stiffness  will  be  unrealistically  lowl 


and  peak  strains  will  be  much  greater  than  they  should  be.  The  samples 
are  vacuum  saturated  with  deaired  sea  water  containing  dissolved  car¬ 
bonates.  Blouin  and  Timian  (1986)  showed  that  the  absence  of  dissolved 
carbonates  in  the  pore  water  results  in  significant  reduction  of  cemen¬ 
tation  strength  in  cemented  samples.  Once  vacuum  saturation  is  complete, 
the  samples  are  allowed  to  come  to  equilibrium  under  about  200  psi  to 
300  psi  pore  pressure.  This  assures  that  any  remaining  air  in  the  samples 
is  driven  into  solution  in  the  pore  water.  An  initial  effective  stress 
approximately  equal  to  the  pore  pressure  is  also  applied.  This  com¬ 
bination  of  initial  pore  pressures  and  effective  stresses  closely  simu 
lates  the  in  situ  condition  of  most  of  the  core  samples  we  tested. 

2.  Undrained  Load-Unload:  The  samples  are  loaded  and  unloaded  statically  in 
an  undrained  condition.  A  peak  axial  stress  on  the  order  of  1  kbar  is 
applied  which  insures  sufficient  strain  to  break  down  the  cementation  in 
all  cemented  materials.  A  correction  must  be  applied  to  the  axial  defor¬ 
mation  to  account  for  the  slight  deformation  in  the  piston  seals,  the  rigid 
vessel,  and  the  pore  pressure  gage  and  line.  This  correction  was  deve¬ 
loped  by  measuring  the  total  axial  displacement  for  a  sample  of  pure  water 
as  a  function  of  pressure  and  then  subtracting  the  portion  of  the  displa¬ 
cement  computed  to  be  due  to  compression  of  the  water.  This  excess 
displacement  due  to  elastic  deformations  in  the  testing  device  is  then 
subtracted  form  the  axial  displacement  data  from  each  test.  This  insures 
that  the  test  specimens  do  not  appear  to  be  softer  than  they  would  be  in 
their  in  situ  condition. 

3.  Consolidation;  Once  each ‘specimen  is  completely  unloaded,  the  pore 
pressure  line  is  opened  and  the  sample  is  reloaded  and  allowed  to  gra¬ 
dually  return  to  its  initial  effective  stress.  This  consolidation  process 
approximates  the  in  situ  process  where  the  material  is  assumed  to  have 
been  liquefied  by  the  dynamic  loading  and  was  then  reconsolidated  under 
its  original  effective  overburden  stress. 

Representative  shock  consolidation  test  data  are  presented  in  Figures  6.2 

through  6.6.  Figures  6.2  and  6.3  show  data  from  saturated  medium  strength  and 
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weak  porous  limestone  samples,  respectively,  and  Figures  6.4  through  6.6  show 
data  from  saturated  uncemented  soil  samples.  A  series  of  three  plots  is 
shown  for  each  test; 

•  total  stress  as  a  function  of  axial  (volume)  strain,  obtained  by 
dividing  the  measured  total  loads  by  the  sample  area; 

•  pore  pressure  as  a  function  of  axial  (volume)  strain,  obtained  from 
the  pore  pressure  gage;  and 

•  effective  stress  as  a  function  of  axial  (volume)  strain;  obtained 

by  subtracting  the  measured  pore  pressure  form  the  total  stress  measure¬ 
ments  at  each  strain  increment. 

We  can  use  Figure  6.2,  showing  data  from  the  low  strength  limestone,  to 
illustrate  the  various  steps  in  the  testing  procedure  and  the  test  results. 
Point  "a"  at  the  start  of  the  undrained  loading  shows  the  initial  conditions, 
i.e.  effective  stress  of  about  300  psi,  pore  pressure  of  about  300  psi  and 

total  stress  of  600  psi.  The  sample  is  loaded  undrained  to  point  "c".  At 

point  "b",  early  in  the  undrained  loading,  the  cementation  in  the  skeleton 
breaks  down.  This  is  clearly  shown  by  the  effective  stress  response  of  Figure 
6.2c.  Prior  to  the  breakdown  of  cementation,  more  than  half  of  the  total 
applied  stress  is  carried  by  the  limestone  skeleton;  but  following  the  cemen¬ 
tation  breakdown  only  5%  of  the  additional  applied  stress  is  carried  by  the 
skeleton.  The  balance  of  the  total  applied  stress  is  carried  by  the  pore 
water. 

From  point  "c",  the  undrained  sample  is  unloaded  to  a  state  of  zero 

stress  at  "e".  As  shown  in  Figure  6.2c,  the  skeleton  is  strongly  hysteretic 

and  unloads  very  rapidly  until  at  "d"  the  effective  stress  drops  to  zero  and  a 
state  of  liquefaction  is  achieved.  Note  that  at  point  "dM  the  overall  volume 
strain  is  still  1.4%  and  the  total  stress  is  still  over  11,000  psi.  From 
point  "d”  to  point  "e"  the  sample  is  liquefied  and  the  total  stress  equals  the 
pore  pressure.  Also  note  that  the  hysteresis  in  the  total  stress  curve  is 
positive,  i.e.  energy  is  being  dissipated,  while  the  apparent  hysteresis  in 
the  pore  pressure  load-unload  cycle  is  negative  as  explained  in  Section  2. 
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At  point  "e"  the  pore  pressure  line  is  opened  and  shortly  thereafter 
total  stress  is  again  applied  and  the  sample  is  allowed  to  drain  and  con¬ 
solidate  until  the  effective  stress  reaches  its  pretest  value  at  point  "f". 

The  strain  during  consolidation  to  the  original  effective  stress  is  nearly  80% 
of  the  strain  reached  during  the  undrained  loading.  The  residual  strain 
following  consolidation,  normalized  by  the  peak  undrained  strain,  is  plotted  as 
a  function  of  dry  density  in  Figure  6.7.  There  is  a  trend  in  the  data  toward 
less  residual  strain  with  increasing  dry  density.  However,  for  samples  with 
dry  densities  of  less  than  100  Ib/ft^,  corresponding  to  porosities  greater 
than  40%  which  are  typical  of  the  Enewetak  materials,  the  average  post  loading 
subsidence  is  85%  of  the  peak  strain  reached  in  the  undrained  loadings.  By 
computing  the  peak  dynamic  strain  fields  form  the  KOA  and  OAK  detonations, 
this  average  consolidation  can  then  be  used  to  estimate  the  portions  of  the 
total  crater  volumes  which  might  be  attributed  to  post  liquefaction  con¬ 
solidation. 

The  other  data  plots  for  the  saturated  limestone  and  soils  of  Figures  6.3 
through  6.6  are  similar  in  character  to  Figure  6.2,  except  that  the  soil  ske¬ 
letons  carry  much  less  of  the  total  stress  and  don't  show  the  early  time 
breakdown  of  the  cementation.  All  the  shock  consolidation  samples,  whether 
cemented  or  uncemented,  showed  very  consistent  and  similar  liquefaction  beha¬ 
vior  early  in  the  unloading  cycle.  We  feel  this  set  of  data  is  a  good  valida¬ 
tion  of  the  theoretical  model  developed  in  Sections  2  and  3. 

6.3  COMPARISON  OF  COMPUTED  AND  MEASURED  UNDRAINEO  UNIAXIAL  RESPONSE. 

The  computational  algorithm  NKOCP,  described  in  Section  3.6,  was  used  to 
compute  the  uniaxial  undrained  response  for  a  saturated  limestone  similar  to 
that  loaded  in  test  A2A6  (see  Table  6.1).  The  comparison  between  the  computed 
response  and  the  test  data  serves  as  a  validation  of  the  numerical  techniques 
and,  as  described  in  Sections  2  and  3,  gives  us  new  insight  into,  and 
understanding  of,  the  processes  controlling  the  two  phase  response. 

The  inputs  to  NKOCP  for  the  saturated  limestone  are  listed  in  Table  6.2 
and  shown  in  Figure  6.8. 
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•  The  constitutive  model  for  the  uniaxial  response  of  the  porous 
skeleton  consists  of  a  tri linear  fit  to  drained  uniaxial  data 
for  porous  Enewetak  limestone  having  nearly  the  same  porosity 

as  sample  A2A6.  This  data  fit  is  shown  in  Figure  6.8  and  consists 
of  an  initial  stiff  loading,  until  the  cementation  breaks  at  a 
strength  of  about  3,150  psi  and  a  strain  of  0.17%,  followed  by 
a  much  softer  loading  to  strains  of  about  5%.  Unloading  follows 
the  same  slope  as  the  initial  cemented  loading. 

•  The  constitutive  model  for  the  sea  water  used  to  saturate  the 
sample  is  that  described  in  Section  3.3  which  is  incorporated  into 
NKOCP  as  option  NW  =  1. 

•  The  constitutive  model  for  the  solid  grains  is  that  described  in 
Section  3.4  using  initial  constants  listed  in  Table  6,2.  These 
constants  were  computed  from  measured  grain  properties  for  the 
Enewetak  materials  described  by  Kim,  Blouin  and  Timian  (1986). 

This  model  was  loaded  using  NKOCP  to  a  peak  total  strain  of  1.587%,  which 
corresponds  to  the  peak  strain  in  the  experiment. 

The  comparison  between  the  calculated  response  and  the  A2A6  test  data  is 
shown  in  Figure  6.9,  parts  a  through  c.  The  good  agreement  between  the  calcu¬ 
lated  response  and  test  data  is  gratifying,  particularly  considering  that  the 
calculated  response  was  developed  from  constitutive  models  for  each  of  the  com¬ 
ponents  which  were  developed  independently  of  the  test  data.  Important 
features  of  the  comparisons  include: 

•  strong  positive  hysteresis  in  the  total  stress  volume  change  cycle 

of  Figure  6.9a,  indicating  that  a  substantial  amount  of  energy  is  dissi¬ 
pated; 

•  apparent  negative  hysteresis  in  both  the  computed  and  measured  pore 
pressure-volume  strain  cycles  of  Figure  6.9b.  As  described  in  the 
discussion  of  Section  2,  this  apparent  negative  hysteresis  is  a  result 
of  rapid  expansion  of  the  solid  grains  as  the  effective  stress  drops  to 
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zero  early  in  the  unloading.  The  expansion  of  the  grains  compresses 
the  pore  water  resulting  in  the  apparent  negative  hysteresis; 

•  strongly  hysteretic  nature  of  the  skeleton  as  evidenced  by  the  effec¬ 
tive  stress-volume  change  comparison  of  Figure  6.9c; 

♦  the  unloading  slope  of  the  skeleton  prior  to  liquefaction  (Figure 
6.9c)  is  significantly  steeper  in  the  calculation  than  in  the  data. 

The  reason  for  this  difference  is  the  assumption  in  the  numerical  model 
that  the  sample  liquefies  as  soon  as  the  radial  effective  stresses  drop 
to  zero.  In  Section  2,  once  the  radial  effective  stresses  reach  zero 
the  skeleton  was  assumed  to  collapse  with  an  immediate  loss  in  axial 
effective  stress.  While  this  may  be  a  valid  model  for  granular  soil, 
it  does  not  appear  to  be  valid  for  the  porous  limestone.  An  apparently 
better  model  would  allow  the  axial  effective  stress  to  return  to  zero 
along  a  uniaxial  unload  path  rather  than  drop  precipitously  once  zero 
radial  effective  stress  is  reached.  Additional  investigation  of  details 
in  the  liquefaction  mechanism  and  refinements  to  the  model  will  be  made 
during  the  coming  year. 

Figure  6.10  shows  the  calculated  total  and  effective  stress  and  pore 
pressure  responses  for  the  limestone  combined  in  a  single  plot.  Note  that  the 
relatively  high  pore  pressures  and  low  initial  Poisson's  ratio  of  the  skeleton 
keep  the  radial  effective  stresses  to  a  minimum.  The  low  peak  radial  effec¬ 
tive  stress  unloads  very  rapidly,  according  to  the  model  assumptions  discussed 
above,  and  causes  liquefaction  very  early  in  the  unloading  cycle. 

The  final  set  of  comparisons  between  calculated  undrained  uniaxial 
response  and  measured  response  is  that  shown  in  Figure  6.11  between  test  A4B6 
and  a  saturated  sand  having  constituent  properties  listed  in  Table  6.3.  The 
skeleton  properties  are  a  fit  to  drained  test  data  for  a  silty  sand  sample 
having  a  porosity  close  to  that  of  sample  A4B6.  Sea  water  and  grain  proper¬ 
ties  are  the  same  as  those  used  in  the  limestone  calculations.  The  com¬ 
parisons  between  the  test  data  and  calculations  are  excellent.  There  are  two 
major  diferences  between  the  soil  response  and  the  limestone  response 
discussed  above: 
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•  The  absence  of  cementation  in  the  soil  sample  results  in  much  less 
of  the  total  load  being  borne  by  the  skeleton; 

•  Even  though  the  skeleton  is  strongly  hysteretie,  the  fact  that  the 
effective  stresss  are  much  less  than  in  the  limestone  means  that 
overall  energy  dissipation  and  hysteresis  in  the  soil  sample  is  much 
less  than  in  the  limestone.  Thus,  the  hysteresis  in  the  soil  samples 
is  minimal  on  the  scale  of  the  total  stresses  and  pore  pressures, 
making  accurate  measurements  of  total  stress  and  pore  pressure 
hysteresis  very  difficult  for  all  the  soil  samples. 
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Table  6.1.  Shock  consolidation  sample  and  test  summary. 


TEST  ID 

MATERIAL 

HOLE  ID/DEPTH 
(ft) 

DRY  DENSITY 
(lb/ft*) 

POROSITY 

(%) 

er/ep* 

M31A6 

Limestone 

KAM2A/245.2 

111.15 

.359 

.65 

M31B6 

Limestone 

KAM2A/362.2 

100.65 

.42 

.76 

M28C6 

Limestone 

KAM2A/432. 1 

92.16 

.469 

.86 

M27B6 

Limestone 

KAM2A/391.3 

115.21 

.336 

.68 

M28B6 

Limestone 

KAM2A/271.Q 

140.91 

.188 

.30 

A2A6 

Limestone 

KAM2A/353.8 

110.85 

.361 

.56 

A4A6 

S/S/G 

0AM1  /  86,3 

95.11 

.452 

.86 

A4B6 

S/S/G 

0AM2  /2S9.8 

103.05 

.406 

.84 

N27A6 

S/S/G 

KAM2  /  59.0 

95.31 

.451 

.79 

N26B6 

S/S/G 

KAM2  /  33.5 

103.22 

.405 

1.0 

A2B6 

S/S/G 

0AM1  /284.9 

100.61 

.42 

.80 

M26A6 

S/S/G 

KAM2  /  11.5 

89.41 

.485 

1.0 

A3A6 

S/S/G 

0AR2  /827.4 

95.14 

.452 

.88 

A386 

S/S/G 

0AR2  /702.4 

104.53 

.398 

.95 

A3C6 

S/S/G 

0AR2  /472.8 

100.18 

.423 

.83 

K17A6 

Beach  Sand 

Beach  Sand 

100.79 

.419 

.73 

M17C6 

Beach  Sand 

Beach  Sand 

103.9 

.40 

.68 

*er/ep  is  the  ratio  of  residual  strain,  er,  following  drained  consolidation, 
to  the  peak  strain,  ep,  attained  during  the  initial  undrained  loading. 
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Table  6.2.  Material  properties  used  for  the  computation  of  undrained  uniaxial 
strain  load-unload  response  of  saturated  porous  Enewetak  limestone. 

Drained  Skeleton  Properties: 

Constrained  loading  modulus  -  Ms]  *  1.86  x  10®  psi  for  ea'  <  .0017 

MS]  »  64,000  psi  for  .0017  4  ea  =  .05 

Constrained  unloading  moduls  -  Msu  =  1.86  x  10®  psi 
K0  Values: 


Loading  - 

ii 

o 

0.02 

V  =  .0196 

€a'  <  0.0017 

n 

o 

0.74 

V  a  .425 

ea*  >  0.0017 

Unloading  - 

K0  = 

0.74 

y  =  .425 

Initial  Porosity:  n  ■  0.361 


Pore  Water  Property: 

Used  the  properties  of  sea  water  described  in  Section  3.3. 

Drain  Properties: 

Used  the  solid  grain  model  described  in  Section  3.4  with  the  following 
properties: 

Initial  grain  density  -  pc  •  0.0002614  lb-s2/in*  (2.78  gr/c»3) 

Initial  wavespeed  -  CQ  »  12,658  ft/s 

Poisson's  ratio  at  low  stress  -  u0  »  0.27 

Slope  of  propagation  to 

particle  velocity  -  S  «  l.S 

Threshold  pressure  -  Pj,  o  72,500  psi 
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Table  6.3.  Material  properties  used  for  the  computation  of  undrained  uniaxial 
strain  load-unload  response  of  saturated  porous  sand. 


Skeleton  Properties: 

Initial  porosity  - 
Constrained  loading  modulus  - 
Constrained  unloading  modulus  - 
Poisson's  ratio  - 


n  a  0.406 
Mgi  =  14,200  psi 
Mgu  =  265,000  psi 
V  =  0.32 


Pore  Water  Properties: 

Used  the  properties  of  sea  water  described  in  Section  3.3. 


Grain  Properties: 

Used  the  solid  grain  model  described  in  Section  3.4  with  the  following 
properties: 

Initial  grain  density  -  p0  •  0.0002614  lb-s2/in<  (2.78  gr/c»2) 

Initial  wavespeed  -  C0  ■  12,668  ft/s 

Poisson’s  ratio  at  low  stress  -  v0  «  0.27 

Slope  of  propagation  to 

particle  velocity  -  S  *  1.5 

Threshold  pressure  -  Pjj  «  72,600  psi 
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Figure  6.1.  Schematic  section  view  of  rigid  high  pressure  oedometer 
used  for  shock  consolidation  tests. 
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ock  Consolidation  (M31B6)  KAM-2A 


consolidation  test  M31B6  on  low  strength  porous  limestone 
stress  response. 


Shock  Consolidation  (M31B6)  KAM-2A 
Depth  -  362.2  ft 


Figure  6.2b.  Shock  consolidation  test  M31B6  on  low  strength  porous  limestone; 
pore  pressure  response. 
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Figure  6.2c.  Shock  consolidation  test  M31B6  on  low  strength  porous  limestone 
effective  strength  response. 
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Figure  6.3a.  Shock  consolidation  test  A2A6  on  medium  strength  porous  limestone 
total  stress  response. 


Shock  Consolidation  (A2A6)  KAM-2A 
Depth  =  353.8  ft 
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Figure  6.3b.  Shock  consolidation  test  A2A6  on  medium  strength  porous  limestone; 
pore  pressure  response. 


Figure  6.3c.  Shock  consolidation  test  A2A6  on  medium  strength  porous  limestone 
effective  stress  response. 


hock  Consolidation  (A4B6)  OAM-2 
Depth  =  259.8  ft 
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Volume  Strain 

Figure  6.4a.  Shock  consolidation  test  A436  on  uncemented  silt,  sand,  gravel 
total  stress  response. 


Shock  Consolidation  (A4B6)  OAM-2 
Depth  =  259. B  ft 


Figure  6.4b.  Shock  consolidation  test  A4B6  on  uncontented  silt,  sand,  gravel 
pore  pressure  response. 


Shock  Consolidation  (A4B6)  OAM 
Depth  =  259.8  ft 


6 
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Figure  6.4c.  Shock  consolidation  test  A4B6  on  uncemented  silt,  sand,  gravel 
effective  stress  response. 


Shock  Consolidation  (A2B6 
Depth  ==  284.9  ft 


Figure  6.5a.  Shock  consolidation  test  A2B6  on  uncemented  silt,  sand,  gravel 
total  stress  response. 


Figure  6.5b.  Shock  consolidation  test  A2B6  on  uncemented  silt,  sand,  gravel 
pore  pressure  response. 


Shock  consolidation  test  A2B6  on  uncemented  silt,  sand,  gravel; 
effective  stress  response. 
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consolidation  test  M27A6  on  uncetnented  silt,  sand,  gravel 
stress  response. 


Figure  6.6b.  Shock  consolidation  test  M27A6  on  uncemented  silt,  sand,  gravel; 
pore  pressure  response. 
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Figure  6.6c.  Shock  consolidation  test  M27A6  on  uncemented  silt,  sand,  gravely 
effective  stress  revxmse. 


residual  volume  strain  following  consolidation  to  peak 
ring  undrained  loading  from  shock  consolidation  tests. 


0.0017  0.05 

Volume  Strain 


Figure  6.8.  Effective  uniaxial  stress-strain  model  for 

calculation  of  the  response  of  saturated  medium 
strength  Enewetak  limestone. 


Consolidation  (A2A6)  KAM-2A 
Depth  =  353.8  ft 


:igure  6.9a.  Comparison  of  total  stress  response  of  saturated  porous  limestone 
with  numerical  simulation. 


Figure  6.9b.  Comparison  of  pore  pressure  response  of  saturated  porous  limestone 
with  numerical  simulation. 


Figure  6.9c.  Comparison  of  effective  axial  stress  response  of  saturated 
porous  limestone  with  numerical  simulation. 


y 


Undrained  uniaxial  strain  load-unload  response  of  saturated 
porous  limestone  computed  using  NKOCP. 
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Figure  6.11a.  Comparison  of  total  stress  response  of  saturated  uncemented 
silt,  sand,  gravel  with  numerical  simulation. 


Shock  Consolidation  (A4B6)  OAM-2 
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Figure  6.11b.  Comparison  of  pore  pressure  response  of  saturated  uncemented 
silt,  sand,  gravel  with  numerical  simulation. 
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Figure  6.11c.  Comparison  of  effective  axial  stress  response  of  saturated 
uncemented  silt,  sand,  gravel  with  numerical  simulation. 


SECTION  7 


GRAIN  DAMAGE  AS  A  FUNCTION  OF  STRESS  PATH  AND  MAGNITUDE 

7.1  INTRODUCTION 

In  Section  2  it  was  shown  that  hysteresis  in  the  material  skeleton 
results  in  liquefaction  following  both  hydrostatic  and  uniaxial  undrained 
loadings  of  porous  saturated  materials.  The  skeleton  hysteresis  is  actually  a 
measure  of  the  non  recoverable  plastic  work  done  on  the  material  during  the 
load  unload  cycle.  There  are  two  principal  plastic  work  mechanisms  acting 
during  the  cycle:  frictional  energy  dissipation  and  grain  breakage. 

Frictional  energy  dissipation  causes  heating  of  the  sample  during  the  load 
unload  cycle,  but  leaves  little  in  the  way  of  post-test  evidence.  Grain 
breakage,  however,  is  easily  recognizable  and  can  be  documented  post-test. 
Blouin,  Martin  and  McIntosh  (1984),  observed  extensive  grain  crushing  during 
uniaxial  and  hydrostatic  drained  loadings  of  Enewetak  beach  sand,  a  well 
rounded  uniform  carbonate  sand. 

In  this  section  we  present  test  results  from  a  systematic  study  of  grain 
crushing  in  a  variety  of  test  types  and  stress  paths  to  several  different  peak 
stresses.  By  relating  the  measured  grain  crushing  to  the  stress-strain  data 
from  the  various  tests,  it  is  hoped  that  we  can  gain  insite  into  the  similari¬ 
ties  and  differences  between  the  different  types  of  loadings  on  the  microscopic 
levels,  and  particularly  into  the  differences  and  similarities  between  the 
microscopic  response  of  granular  materials  under  shear  and  compressive 
loadings.  The  data  presented  in  this  section  will  support  our  investigation 
and  analysis  of  grain  crushing  during  the  coming  year. 

7.2  TEST  DESCRIPTIONS 

A  series  of  nine  drained  compression  tests,  listed  in  Table  7.1,  was  run 
on  well  rounded  Enewetak  beach  sand  samples  which  had  been  sieved  and  sorted 
so  that  all  the  pretest  material  was  retained  on  the  number  40  sieve  and  had  a 
grain  size  ranging  from  0.425  to  0.60  mm.  All  test  samples  were  cylindrical 
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with  a  nominal  diameter  of  2  in  and  length  of  4  in.  All  samples  were  loaded 
in  a  triaxial  vessel  along  one  of  four  stress  paths: 

1.  hydrostatic  compression  (HY)  to  pressures  of  1,000  psi  and 
10,000  psi; 

2.  uniaxial  strain  compression  (K0)  to  mean  stresses  of  approximately 
1,000  psi  and  10,000  psi; 

3.  triaxial  compression  (TXC)  to  mean  stresses  of  approximately 
1,000  psi,  10,000  psi  and  15,300  psi;  and 

4.  triaxial  shear  at  constant  mean  stresses  (TXS)  of  approximately  1,000 
and  10,000  psi. 

The  stress  paths  are  shown  on  a  plot  of  stress  difference  vs.  mean  stress 
in  Figures  7.1a  and  b.  The  hydrostatic  loading  (HY)  simply  moves  along  the 
mean  stress  axis  with  no  shear  stress  generated.  The  K0  uniaxial  loading 
follows  a  stress  path  governed  by  the  sample  as  it  deforms  axially  with  no 
lateral  strain  permitted.  The  triaxial  compression  test  (TXC)  follows  a 
prescribed  stress  path  as  the  sample  is  loaded  axially  under  a  constant  con¬ 
fining  pressure.  The  confining  pressure  (either  630  psi  or  8,400  psi)  was 
computed  so  that  the  final  stress  state  would  approximately  equal  that  from 
the  corresponding  K0  test.  Two  TXC  tests  were  run  at  8,400  psi  confining 
pressure;  test  Q60E7  to  a  mean  stress  of  about  10,000  psi  and  test  U23AB7  to  a 
mean  stress  of  15,350  psi.  The  triaxial  shear  tests  (TXS)  were  run  at 
constant  mean  pressures  to  approximately  the  same  final  stress  states  as  the 
corresponding  Kq  and  TXC  tests. 


7.3  TEST  RESULTS 

Post-test  microscopic  views  of  representative  grain  samples  showing  the 
grain  damage  are  presented  in  Figures  *fc2a  through  c.  In  Figure  7.2a,  a  pre¬ 
test  sample  is  compared  to  a  post-test  sample  from  test  U23AB7,  the  TXC  test 
to  15,350  psi  mean  stress.  The  uniform  well  rounded  nature  of  the  grains  is 
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apparent  in  the  virgin  material.  Also  note  the  presence  of  microporosity 
visible  in  some  of  the  individual  grains.  Extensive  grain  fracturing  and 
crushing  is  evident  in  the  post-test  sample,  with  only  a  few  grains  apparently 
having  survived  intact.  Most  of  the  post-test  material  has  been  broken  into  a 
wide  variety  of  grain  sizes. 

Figure  7.2b  shows  post- test  microscopic  views  of  material  from  the  four 
tests  run  to  a  peak  mean  stress  of  about  1,000  psi.  Qualitatively,  there 
appears  to  be  little  difference  between  the  four  samples^  In  all  cases  only 
minor  grain  damage  is  apparent,  with  a  modest  number  of  grains  broken  into 
large  pieces  and  creation  of  only  a  small  amount  of  finer  grain  sizes. 

Microscopic  views  of  the  material  from  the  four  10,000  psi  mean  stress 
tests  are  compared  in  Figure  7.2c.  Grain  damage  is  substantially  more  per¬ 
vasive  than  that  from  the  1,000  psi  tests,  but  again  there  does  not  appear  to 
be  a  significant  difference  in  damage  between  the  four  tests.  The  10,000  psi 
results  show  a  more  well  graded  post-test  grain  distribution  than  the  1,000 
psi  tests,  with  numerous  medium  sized  and  small  sized  grains  in  evidence. 

The  post-test  grain  size  distributions  from  the  nino  tests  are  plotted  in 
parts  "aM  of  Figures  7.3  through  7.11,  and  listed  in  Table  7.1.  The  distribu¬ 
tions  are  plotted  as  percent  by  weight  within  a  given  size  interval  as  a  func¬ 
tion  of  the  log  of  the  interval.  Thus  each  interval  is  represented  by  an 
equal  width  and  the  area  within  each  interval  is  proportional  to  the  weight  of 
material  having  grain  sizes  within  the  interval.  Gradations  were  not  run  on 
material  passing  the  number  200  sieve  (smaller  than  0.076  mm).  For  con¬ 
venience  all  such  material  was  lumped  into  a  single  grain  interval  equal  in 
width  to  the  others,  but  delineated  by  a  dashed  line  such  as  shown  in  Figure 
7.lla, 


Figures  7.12  and  7.13  compare  the  post-test  grain  size  distributions  from 
the  1,000  psi  and  10,000  psi  mean  stress  tests,  respectively.  The  samples 
loaded  to  1,000  psi  mean  effective  stress  sustained  only  minor  grain  damage, 
with  between  3%  and  7%  of  the  pretest  grains  being  broken  into  a  smaller  size 
interval.  Most  of  the  broken  grains  only  moved  to  the  next  smaller  interval, 
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suggesting  that  the  original  grains  sustained  only  one  or  two  fractures.  This 
conclusion  is  supported  qualitatively  by  the  post-test  microscopic  views  of 
Figure  7.2b.  The  material  loaded  in  the  hydrostatic  test  received  the  least 
amount  of  damage,  with  only  2.8%  of  the  grains  broken  to  a  smaller  size.  The 
three  samples  loaded  in  combined  hydrostatic  and  shear  stress  sustained 
somewhat  more  damage  (between  5.3%  and  6.7%  broken  to  a  smaller  size).  The 
amount  of  damage  to  these  three  samples  is  very  similar  in  all  cases  and 
appears  to  be  independent  of  the  stress  path  used  to  reach  the  final  state  of 
stress. 

The  comparison  of  grain  damage  for  the  10,000  psi  mean  stress  tests  shown 
in  Figure  7.13  shows  a  moderate  to  high  degree  of  damage,  36%  to  49%  of  the 
grains  broken  into  a  smaller  size  interval.  In  all  cases,  grain  breakage  is 
much  more  extensive  than  breakage  from  the  1,000  psi  loadings,  with  a  well 
graded  distribution  of  broken  grains  being  created.  Again,  the  material 
loaded  in  hydrostatic  compression  sustained  the  least  damage,  with  36.4%  of 
the  grains  reduced  in  size.  Damage  form  the  combined  hydrostatic  and  shear 
loadings  ranged  from  42.2%  for  the  TXC  test  to  49.1%  for  the  K0  test.  The 
differences  between  these  three  is  not  thought  to  be  significant,  though  addi¬ 
tional  testing  would  be  required  to  substantiate  this. 

While  the  material  loaded  in  hydrostatic  compression  sustained  loss 
damage  in  both  test  series,  the  difference  in  damage  between  the  hydrostatic 
tests  and  the  shear  tests  is  quite  small.  Additional  tests  at  10,000  psi  mean 
stress  but  with  higher  shear  stress  components  are  planned  to  better  define 
the  influence  of  shear  stress  on  grain  damage.  Severe  grain  damage  was 
sustained  by  the  triaxial  sample  loaded  under  the  same  8,400  psi  confining 
pressure  a?  the  10,000  psi  mean  stress  sample,  but  carried  to  a  maximum  mean 
stress  of  15,350  psi  and  a  maximum  stress  difference  of  20,900  psi.  Nearly 
80%  of  the  grains  were  reduced  in  size,  with  over  half  the  total  post-test 
weight  having  grain  sizes  of  less  than  .106  mm. 

In  order  to  briefly,  but  quantitatively,  describe  the  grain  damage  we 
developed  a  technique  which  uses  the  pretest  and  post-test  log  mean  grain 
sizes  to  compute  a  grain  damage  factor,  Df,  given  by 
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where 


D0  *  pre-test  log  mean  grain  size 

Dp  =  post-test  log  mean  grain  size 

Dra  »  minimum  measured  grain  size  interval 


The  log  mean  grain  size  represents  the  center  of  mass  of  the  grain  size  and 
weight  fraction  semi-log  histograms  in  terms  of  grain  size.  Thus  all  of  the 
pretest  samples  have  a  mean  grain  size  of  0.50  mm,  while  the  post-test  mean 
grain  sizes  range  from  0.495  mm  for  the  1,000  psi  HY  test  (Y4A7)  to  0.177  mro 
for  the  15,350  psi  mean  stress  TXC  test  (U23A87). 


The  grain  damage  factor  is  a  measure  of  both  the  amount  and  severity  of 
the  grain  damage.  A  sample  with  most  of  its  grains  crushed  to  a  fine  powder 
will  have  a  much  larger  grain  damage  factor  than  a  sample  with  most  of  its 
grains  broken  into  only  slightly  smaller  pieces.  As  shown  in  Table  7.1,  the 
grain  damage  factor  ranges  from  0.011  for  the  1,000  psi  HY  sample  to  0.743  for 
the  15,350  mean  stress  TXC  sample.  The  maximum  grain  damage  factor  is  1.0. 


ihe  stress-strain  response  of  the  material  in  the  nine  loading  tests  is 
presented  in  parts  b  through  h  of  Figures  7.3  through  7.11.  These  test 
results  will  be  used  in  future  analysis  to  compare  the  plastic  work  done  on 
the  samples  to  grain  damage,  and  to  evaluate  grain  damage  as  a  function  of 
other  response  parameters. 
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Stress  Difference  07-  cr  (psi) 


Figure  7.1b.  Stress  paths  to  approximate  mean  stress  of  10,000  psi. 


168 


a.  Virgin  sand 


b.  Post  test  sand 


0  Scale  (mm)  5 

I - 1 - 1 - 1 - T - 1 

Figure  7.2a.  Microscopic  comparison  of  virgin  beach  sand  with 
post-test  sample  loaded  in  triaxial  compression  to 
a  mean  stress  of  15,300  psi. 


c.  Triaxial  compression 


d.  Triaxial  shear  at  constant  mean  stres: 


Scale  (mm) 


Figure  7.2b.  Microscopic  comparison  uf  beach  sand  sample.,  loaded  along  various 
stress  paths  to  a  mean  stress  of  approximately  1000  psi. 


190 


c.  Tr i axial  compression  <j.  Triaxial  shear  at  constant 

mean  stress 

0  Scale  (mm)  5 

I - T - T - ! - 1 - 1 

Figure  7.2c.  Microscopic  comparison  of  beach  sand  samples  loaded  along  various  stress 
paths  to  a  mean  stress  of  approximately  10,000  psi. 
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Figure  7.3a.  Post-test  grain  size  and  wei 
hydrostatic  drained  loading 
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Figure  7.3b.  Pressure-volume  strain  data,  hydrostatic  loading  to  1,000  psi. 
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Figure  7.4a.  Post  test  grain  size  and  weight-fraction  distribution  following  K( 
drained  loading  to  approximately  1000  psi  mean  stress. 


Axial  stress-axial  strain  data,  Kq  loading  to  approximately  1000  psi 
mean  stress. 
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Figure  7.4c.  Axial  stress-confining  stress  data,  K0  loading  to  approximately  1000 
psi  mean  stress. 


Volume  Strain  {%) 

Mean  stress-volume  strain  data,  K*,  loading  to  approximately  1000  psi 
mean  stress. 


Figure  7.4e.  Stress  difference-mean  stress  data,  Kc  loading  to  approximately  1000 
psi  mean  stress. 
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Stress  difference-volume  strain  data,  K0  loading  to  approximately  1000 
psi  mean  stress. 
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Figure  7.5a.  Post-test  grain  size  and  weight  fraction  distribution  following 
triaxial  compression  to  approximately  1000  psi  mean  stress. 


Hydrostat  (G5B7)  to  630  psi 
Beach  Sand 


Volume  Strain  (%) 

Figure  7.5b.  Pressure-volume  strain  data  to  630  psi  prior  to  triaxial  compression 
approximately  1G00  psi  mean  stress. 


5C7)  from  630  psi 
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ess  difference-axial  and  radial  strain  data,  triaxial  compression 
approximately  1000  psi  mean  stress. 
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Figure  7.5e.  Stress  difference-volume  strain  data,  triaxial  compression  to  approxi 
mately  1000  psi  mean  stress. 


G5C7)  from  630  psi 
Beach  Sand 


stress-volume  strain  data,  triaxial  compression  to  approximately 
psi  mean  stress. 


Triax  (G5C7)  from  630  psi 
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Figure  7.5g.  Radial  strain-axial  strain  data,  triaxial  compression  to  approximately 
100Q  psi  mean  stress. 
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Figure  7.5h.  Strain  difference-volume  strain  data,  triaxial  compression  to  approxi¬ 
mately  1000  psi  mean  stress. 
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Figure  7.6a.  Post  test  grain  size  and  weight-fraction  distribution  following 
triaxial  shear  at  constant  mean  stress  of  approximately  1000  psi 
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r.6c.  Mean  stress-volume  strain,  hydrostatic  loading  and  TXS  loading  at 
constant  mean  stress  of  approximately  1000  psi. 


Hydrosfcat  w/  Triax  (G6A7)  1000  psi 

Beach  Sand 
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Figure  ?*6e.  Stress  difference- strain  difference  data,  hydrostatic  loading  and  TXS 
loading  at  constant  wean  stress  of  approximately  1000  psi. 


Triax  (G6A7)  1000  psi 
Beach  Sand 


igure  7.6f.  Stress  difference-volume  strain  data,  hydrostatic  loading  and  TXS  load 
ing  at  constant  mean  stress  of  approximately  1000  psi. 
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Figure  7.6g.  Radial  strain-axial  strain  data,  hydrostatic  loading  and  TXS  loading 
at  constant  meant  stress  of  approximately  1000  psi. 


Hydrostat  w/  Triax  (66A7)  1000  psi 
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Figure  7.6h.*  Strain  difference-volume  strain  data,  hydrostatic  loading  and  TXS  load¬ 
ing  at  constant  mean  stress  of  approximately  1000  psi. 
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Figure  7.7b.  Pressure-volume  strain  data,  hydrostatic  loading  to  10,000  psi. 
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Figure  7.8a.  Post  test  grain  size  and  weight-fraction  distribution  following  Kq 
drained  loading  to  18,000  psi  mean  stress. 
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Confining  Stress  (psi) 

Figure  7.8c.  Axial  stress-confining  stress  data,  K0  loading  to  approx 
psi  mean  stress. 
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Figure  7.8d.  Mean  stress-volume  strain  data,  Kq  loading  to  approx 
psi  mean  stress. 
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Figure  7.8e.  Stress  difference-mean  stress  data, 
psi  mean  stress. 
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Grain  Size  (mm) 

Figure  7.9a.  Post  test  grain  size  and  weight-fraction  distribution  following 
triaxial  compression  to  approximately  10,000  psi  mean  stress. 


Figure  7.9b.  Pressure-volume  strain  data  to  8400  psi  prior  to  triaxial  compression 
to  approximately  10,000  psi  mean  stress. 
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approximately  10,000  psi  mean  stress. 


Triax  (G5E7)  from  8400  psi 
Beach  Sand 


Figure  7.9d.  Stress  difference-strain  difference  data,  triaxial  compression  to 
approximately  10,000  psi  mean  stress. 
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Figure  7. 9e.  Stress  diff erence- volume  strain  data,  triaxial  compression  to  approx¬ 
imately  10,000  psi  mean  stress. 


Figure  7.9f.  Kean  stress-volume  strain  data,  triaxial  compression  to  approximately 
10,000  psi  mean  stress. 
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Figure  7.9g.  Radial  strain-axial  strain  data,  triaxial  compression  to  approximately 
10,000  psi  mean  stress. 
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Figure  7.10a.  Post  test  grain  size  and  weight- fraction  distribution  following 

triaxial  shear  at  constant  mean  stress  of  approximately  10,000  psi 
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Figure  7.10b.  Axial  stress-axial  strain  data,  hydrostatic  loading  and  TXS  loading 
at  constant  mean  stress  of  approximately  10,000  psi. 
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Figure  7.10c.  Mean  stress-volume  strain  data,  hydrostatic  loading  and  TXS  loading 
at  constant  mean  stress  of  approximately  10,000  psi. 
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Stress  difference-axial  and  radial  strain  data,  hydrostatic  loading 
and  TXS  loading  at  constant  mean  stress  of  approximately  10,000  psi 
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Figure  7.10e.  Stress  difference-strain  difference  data,  hydrostatic  loading  and  TXS 
loading  at  constant  mean  stress  of  approximately  10,000  psi. 
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rigure  7.10f.  Stress  difference-volume  strain  data,  hydrostatic  loading  and  TXS  load 
ing  at  constant  mean  stress  of  approximately  10,000  psi. 
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Figure  7.10g.  Radial  Strain-axial  strain  data,  hydrostatic  loading  and  TXS  loading 
at  constant  mean  stress  of  approximately  10,000  psi. 
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igure  7.10h.  Strain  difference- volume  strain  data,  hydrostatic  loading  and  TXS 
loading  at  constant  mean  stress  of  approximately  10,000  psi. 
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Figure  7.11a.  Post  test  grain  size  and  weight  fraction  distribution  following 
triaxial  compression  to  15,350  psi  mean  stress. 
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Figure  7.11c.  Stress  difference-axial  and  radial  strain  data,  triaxial  compression 
to  15,350  psi  mean  stress. 
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Figure  7. lid.  Stress  difference-strain  difference  data,  triaxial  compression  to 
15,350  psi  mean  stress. 
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Figure  7*lle.  Stress  difference-volume  strain  data,  triaxial  compression  to  15,350 
psi  mean  stress. 
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Figure  7.ilf.  Mean  stress-volume  strain  data,  triaxial  compression  to  15,350  psi 
mean  stress. 
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Figure  7.11g.  Radial  strain-axial  strain  data,  triaxial  compression  to  15,350  psi 
mean  stress. 


Triax  (U23B7) 

Beach  Sand,  Confining  Pressure  =  8400  psi 


247 


Figure  7.11h.  Strain  difference-volume  strain  data,  triaxial  compression  to  15,350 
psi  mean  stress. 


Size  Distributions 


Figure  7.12.  Comparison  of  post  test  grain  size  and  weight  fraction  distributions 
from  samples  loaded  to  1000  psi  mean  stress. 
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Figure  7.13.  Comparison  of  post  test  grain  size  and  weight  fraction  distribution 
from  samples  loaded  to  10,000  psi  and  15,350  psi  mean  stress. 


SECTION  8 


LABORATORY  INVESTIGATION  OF  FLUID  FLOW  THROUGH  POROUS  MEDIA 

8.1  INTRODUCTION 

A  key  component  of  multiphase  modeling  is  the  equation  of  motion  of  the 
fluid  within  the  porous  skeleton.  Nearly  all  previous  work  in  this  area  is 
based  on  work  performed  by  Biot  (1956,  1962a,  1962b).  Biot  determined  theore¬ 
tical  relationships  for  nonsteady  state  laminar  flow  through  circular  and  flat 
ducts. 

Ward  (1964)  conducted  a  series  of  experimental  flow  tests  and  based  upon 
analysis  of  these  tests,  he  proposed  an  empirical  equation  which  is  applicable 
for  both  laminar  and  turbulent  flow  in  porous  media.  However,  Ward's  work  was 
limited  to  steady  state  flow. 

A  proposed  viscous  friction  equation  combining  our  generalization  of 
Biot’s  theoretical  equation  and  Ward's  empirical  equation  is  given  by: 

04  ■  It  flu  ♦  Hi  (iSu)  +  SI  r  W4 
1  k  1  K*  "  (8.1) 

where  «  viscous  friction  force  between  the  pore  fluid  and  the  soil 

skeleton  per  unit  volume  of  pore  fluid 

Yf  ■  unit  weight  of  the  pore  fluid 

Pf  *  mass  density  of  the  pore  fluid 

n  ■  porosity 

k  »  coefficient  of  permeability 
r  «  mass  increment  factor 

frf  =  Ward's  empirical  constant  which  is  a  function  of  both 
skeleton  and  fluid  properties 
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w-j  =  apparent  fluid  velocity 


w-;  *  apparent  relative  fluid  acceleration 

This  equation  requires  validation  and/or  modification  based  on  laboratory 
experiments.  The  laboratory  investigation  of  fluid  friction  in  nonsteady 
state  flow  and  in  the  turbulent  regime  is  an  important  part  of  the  experimen¬ 
tal  portion  of  our  research.  A  research  plan  was  formulated  which  includes 
the  following  goals: 

1.  Design  and  fabrication  of  a  high  pressure  laboratory  fluid  flow  device. 

2.  Experimental  verification  of  Biot's  theoretical  equations. 

3.  Validation  of  Ward's  flow  model. 

4.  Experimental  tests  to  verify  and/or  modify  Equation  8.1. 

5.  Development  of  relationships  for  nonuniform  flow  paths  under  steady 
state  flow. 

6.  Development  of  relationships  for  nonsteady  state  flow  through  nonuniform 
paths  in  the  turbulent  flow  regime. 

8. 2  RESEARCH  ACCOMPLISHMENTS 

8.2.1  Design  and  Fabrication 

The  design  process  for  the  fluid  flow  device  included  determination  of 
equipment  and  testing  criteria  to  accomplish  the  goals  listed  above.  These 
criteria  include  sample  size  requirements  for  both  uniform  and  nonuniform 
ducts,  pressure  and  pressure  gradient  conditions,  test  control  and  data 
acquisition  requirements. 

It  was  determined  that  the  device  oust  be  capable  of  testing  samples 
2.0  inches  in  diameter  and  9.0  inches  in  length.  Other  capabilities  include 
the  ability  to  pressurize  the  fluid  to  5,000  psi,  the  ability  to  measure  and 
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control  the  rate  of  flow  through  the  sample,  the  ability  to  measure  the 
pressure  and  temperature  on  the  high  and  low  pressure  sides  of  the  sample, 
and  the  ability  to  change  either  the  fluid  or  sample  easily. 

Based  on  the  design  criteria,  a  flow  device  was  designed  and  constructed 
utilizing  components  of  a  5,000  psi  hydraulic  cylinder.  A  schematic  of  the 
device  is  shown  in  Figure  8.1.  A  comp  iter  controlled  servo  .valve  provides 
hydraulic  oil  pressure  to  the  rod  side  of  the  piston  which  pressurizes  the 
large  fluid  reservoir  shown.  The  computer  control  allows  various  loading  con¬ 
ditions  to  be  programmed  including  specified  steady  state  conditions  and  a 
variety  of  nonsteady  state  flow  conditions  including  constant  velocity, 
constant  acceleration,  and  variable  acceleration.  By  monitoring  flow  under 
those  various  conditions  the  terms  of  Equation  8.1  can  be  investigated  inde¬ 
pendently.  Flow  is  measured  by  monitoring  the  position  of  the  loading  piston. 
Based  on  the  position  of  the  piston  at  each  time  step  and  the  area  of  the 
reservoir,  a  simple  calculation  yields  the  flow  through  the  sample. 

Transducers  measure  the  pressure  and  temperature  on  both  the  high  and  low 
pressure  sides  of  the  sample  from  which  the  pressure  and  temperature  gradients 
can  be  calculated. 

8.2.2  Steady  State  Flow  Tests,  Circular  Duct 

Steady  state  flow  tests  were  conducted  through  both  a  circular  and  a  flat 
duct.  The  circular  duct  has  a  diameter  of  0.0625  inches  and  a  length  of  9.125 
inches. 

Oata  output  from  each  of  a  series  of  steady  state  flow  tests  is  in  the 
overall  form  of  displacement  and  pressure  difference  time  histories.  Figures 

8.2  and  8.3  show  typical  output  from  a  circular  duct  test.  A  velocity  value 
is  measured  from  thj  displacement  time  history  by  finding  the  slope  of  the 
line.  The  slope  of  the  line  is  taken  between  two  time  intervals  when  the 
pressure  difference  is  constant. 

Figure  8.4  shows  the  results  of  36  steady  state  flow  tests  using  kerosene 
as  the  pora  f >uid.  The  velocities  range  in  value  from  2.5  inches  pc'*  second 


to  2,350  inches  per  second.  The  pressure  gradient,  which  is  equal  to  the 
overall  pressure  difference  divided  by  the  duct  length,  is  plotted  as  a  func¬ 
tion  of  the  velocity.  Figure  8.5  is  an  enlargement  of  Figure  8.4  below  180 
in/s.  The  initial  portion  of  the  line  defined  by  the  points  below  about  60 
in/s  is  clearly  linear  which  corresponds  to  the  laminar  flow  regime.  Figure 
8.4  shows  that  the  points  transition  from  a  linear  relationship  into  a  nonli¬ 
near  relationship  where  the  pressure  gradient  is  proportional  to  the  square  of 
the  velocity.  The  latter  corresponds  to  the  turbulent  flow  regime  and  the 
transition  region  between  the  linear  and  turbulent  regimes  is  the  transitional 
flow  regime. 

The  governing  equation  of  motion  for  pore  fluid  flow  derived  from  Biot's 
work  (see  previous  references)  expresses  the  pore  pressure  gradient,  ir,f,  as 


»  pf  iii  +  0-j  (8.2) 

where  ir,  \  *  pressure  gradient 

U(  «  acceleration  of  the  pore  fluid 

For  steady  state  flow  through  a  circular  or  flat  duct.  Equation  8.2  can  be 
combined  with  Equation  8.1  to  produce  the  form 

v,{  «  &  *2  (8.3) 

k  k* 


Dividing  both  sides  of  w  produces  an  equation  of  the  form 


where 


y  =  a  ♦  bx 
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Replotting  the  steady  state  data  in  the  form  of  the  pressure  gradient  divided 
by  the  relative  fluid  velocity  vs  the  relative  fluid  velocity  allows  a  linear 
regression  to  be  performed.  The  intercept  can  be  equated  with  "a"  and  the 
slope  is  equal  to  "b". 


Knowing  the  unit  weight  of  the  fluid  allows  the  calculation  of  the  coef¬ 
ficient  of  permeability  by 


k  =  Jf  (8.4) 

Once  the  coefficient  of  permeability  has  been  calculated.  Ward's  constant, 

/3f,  can  then  be  evaluated  by  solving  the  equation. 

fit  *  k*  b  (8.5) 

Figure  8.6  shows  the  steady  state  data  from  the  circular  duct  replotted 
with  the  pressure  gradient  divided  by  the  velocity  plotted  as  a  function  of 
the  velocity.  A  linear  relationship  is  clearly  evident.  A  regression  analy¬ 
sis  provides  an  intercept  of  1.91  x  10“3  Ibs-s/in*  and  a  slope  of  2.069  x  10~8 
lb-s/in8.  Knowing  the  unit  weight  of  kerosene  is  approximately  2.914  x  10”2 
lb/in3  allows  the  calculation  of  the  coefficient  of  permeability  from  Equation 
8.4.  For  the  circular  duct,  the  experimental  value  of  k  is  15.26  inches  per 
second.  Kim  and  Blouin  (1964)  showed  that  k  can  also  be  calculated  for  the 
circular  duct  theoretically  from  the  formula 

k  =  (8.6) 

8|i 

where  a  »  the  radius  of  the  duct 

p  s  dynamic  viscosity  of  the  pore  fluid 
n  =  porosity 

Tf  =  unit  weight  of  the  pore  fluid 
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Subs t i cut i ng  the  appropriate  values  into  Equation  8.6  yields  a  theoretical 
value  of  permeability  equal  to  13.45  in/s.  This  value  is  satisfyingly  close 
to  the  15.26  in/s  measured  experimentally. 

Ward's  constant  can  also  be  calculated  using  Equation  3.5  and  the 
measured  coefficient  of  permeability,  fif  works  out  to  be 
8.082  x  10“5  lb-si-Vin4-5. 

8.2.3  Constant  Acceleration  Flow  Tests,  Circular  Duct 

Figures  8.7  and  8.8  present  the  displacement  and  pressure  difference  time 
histories  for  a  constant  acceleration  flow  test  through  the  circular  duct.  An 
equation  describing  the  displacement  curve  is 

w  =  4.30  t  +  .526  t2  (8.7) 

which  can  be  differentiated  to  produce  the  velocity  and  acceleration  equations. 
These  equations  are 


w  =  4.30  +  1.05  t 


(8.8) 


w  =  1.05 


To  compare  this  constant  acceleration  dynamic  test  data  to  steady  state 
test  data,  we  computed  flow  velocities  and  pressure  gradients  at  40  and  80 
seconds.  As  shown  in  Figure  8.9,  the  dynamic  test  data  agrees  with  the  steady 
state  data,  indicating  that  the  magnitude  of  the  fluid  acceleration  is  too  low 
to  produce  a  measurable  influence  of  fluid  inertia  and  mass  increment  factor 
on  the  fluid  friction.  Attempts  will  be  made  to  increase  accelerations  so  that 
these  influences  can  be  measured. 
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8.2.4  Steady  State  Flow  Tests,  Flat  Duct 


As  a  parallel  effort  to  the  circular  duct,  we  also  conducted  a  set  of 
steady  state  flow  tests  through  a  flat  duct.  The  flat  duct,  as  shown  in 
Figure  8.10a,  has  a  height  of  0.04  inches,  a  width  of  0.28  inches,  and  a 
length  of  8.71  inches.  The  flat  duct  dynanic  flow  equations  developed  by  Biot 
assume  that  the  width  of  the  duct  is  infinite  so  that  there  will  be  no 
influence  on  fluid  friction  due  to  the  side  walls.  Our  flat  duct,  however, 
has  finite  width  with  an  aspect  ratio  of  seven.  To  check  the  accuracy  of  the 
steady  state  test  results,  as  well  as  the  influence  of  the  side  wall  friction, 
we  developed  the  method  of  calculating  the  equivalent  flow  system  in  Appendix 
E.  A  total  of  39  steady  state  flow  tests  through  the  flat  duct  were  conducted. 
These  results  are  plotted  in  Figure  8.11  in  terms  of  flow  velocity  vs  pressure 
gradient,  with  an  expanded  view  of  the  low  velocity  region  of  Figure  8.11 
given  in  Figure  G.12.  As  with  the  circular  duct,  we  have  replotted  the  steady 
state  data  in  terms  of  the  pressure  gradient  divided  by  the  fluid  velocity  vs 
the  fluid  velocity  in  Figure  8.13.  A  least  squares  fit  to  data  yields  the 
coefficient  of  permeability,  k  ■  11.88  in/s,  and  Ward's  constant, 
t if  o  5.0  x  10"6  lb-sl*5/in4’5.  The  theoretical  value  of  k,  computed  in 
Appendix  6,  is  11.92  in/s  which  is  very  close  to  the  experimental  value  of 
11.05  in/s. 

8.2.5  Constant  Acceleration  Flow  Tests,  Flat  Ouct 

Three  sets  of  constant  acceleration  dynamic  flow  tests  through  the  flat 
duct  were  conducted  at  three  different  fluid  accelerations.  Figures  8.14, 

8.15,  and  8,16  show  the  dynamic  test  results  for  fluid  accelerations  of  27.1, 
42.7,  and  63.1  in/s2,  respectively.  For  each  test,  time  histories  of  actual 
ram  displacements,  high  pressures  at  the  duct  inlet  and  low  pressure  at  the 
duct  outlet  are  presented. 

To  evaluate  the  influence  of  fluid  inertia  and  mass  increment  factor  on 
the  fluid  friction,  we  computed  flow  velocities  and  pressure  gradients  at 
various  times  and  compared  them  to  the  steady  state  data  as  shown  in  Figure 
8.17.  The  dynamic  data  presented  in  Figure  8.17  are  slightly  above  the  fit  to 
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the  steady  state  data  in  the  turbulent  flow  regime.  The  deviation  of  the 
dynamic  data  from  the  steady  state  fit  increases  slightly  with  increasing 
fluid  acceleration,  but  higher  fluid  accelerations  must  be  achieved  in  order  tc 
make  confident  measurements  of  the  inertial  and  mass  increment  influences. 
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:igure  8.1.  Schematic  section  view  of  servo  controlled  fluid  flow 
device. 
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Figure  8.2.  Displacement  time  history  for  steady  state  flow  through  circular  duct. 
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Steady  State  Flow  through  Circular  Duct 
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Figure  8.4.  Steady  state  data  showing  laminar,  transition  and  turbulent  regimes 


Figure  8.5.  Expanded 
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Figure  8.7.  Displacement  time  history  from  constant  acceleration  test 


Constant  Acceleration  Flow  through  Circular  Duct 


Figure  8.8.  Pressure  difference  time  history  from  constant  acceleration  test. 


Constant  acceleration  dynamic  flow  test  through  a  circular  duct,  compared 
with  steady  state  flos*  data. 


Figure  8.10.  Cross  section  view  of  flat  duct  and  subdivisions  of  cross  section  used  to 
compute  equivalent  coefficient  of  penneabil ity . 
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Figure  8.11.  Steady  state  data  showing  laminar,  transition  and  turbulent 
regimes,  flat  duct  flow  test. 
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Expanded  view  of  Fi 
fiat  duct  flow  test 
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ady  state  data  replotted  in  order  to  measure  the  coefficient 
permeability  and  Ward's  constant,  flat  duct  flow  test. 


r 


r- 

u 

CD 

C\J 


CJ 


oo 


LO 


CM 


O) 


to 


o 

0) 

CO 


o 


4 

(SJIW)  uipy  | 

"5s 


272 


Figure  8.14a.  Displacement  time  history  from  constant  acceleration  dynamic  flow  test 
through  a  flat  duct,  acceleration  =  27.1  in/sec2. 


1.14b.  High  pressure  time  history  from  constant  acceleration  dynamic  flow  test 
through  a  flat  duct*  acceleration  -  27.1  in/sec2. 
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Figure  8.14c.  Low  pressure  time  history  from  constant  acceleration  dynamic  flow  test 
through  a  flat  duct,  acceleration  =  27.1  in/sec2. 
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Figure  8,15a.  Displacement  time  history  from  constant  acceleration  dynamic  flow  test 
through  a  flat  duct,  acceleration  =  42.7  in/sec2. 
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Figure  8.15b.  High  pressure  time  history  from  constant  acceleration  dynamic  flow  test 
through  a  flat  duct,  acceleration  =42.7  in/sec2. 


15c.  Low  pressure  time  history  from  constant  acceleration  dynamic  fl 
through  a  flat  duct,  acceleration  =  42.7  in/sec2. 


displacement  in  duct,  multiply  the  ram  displacement 


Figure  8.16a.  Displacement  time  history  from  constant  acceleration  dynamic  flow  test 
through  a  flat  duct,  acceleration  =  63.1  in/sec2. 


Figure  8.16b.  High  pressure  time  history  from  constant  acceleration  dynamic  flow 
test  through  a  flat  duct,  acceleration  =63.1  in/sec2. 
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APPENDIX  A 

LISTING  OF  PROGRAM  NKOCP 
C 

C  REVISION  :  JULT  1,  1987  (ALTERED  TO  RUN  ON  IBM-PC-AT) 

C  INCLUDE  INITIAL  STRESSES  AND  UNLOADING 

C  REVISION  :  JULY  5,  1987  (NO  CHANGES  NEEDED  TO  RUN  ON  VAX  II) 

C 

COMMON  /SKTON/  RR(IOO) , CM (100) 

COMMON  /KGDAT/  EKQ (100) ,RKO (100) 

COMMON  /UNLDP/  A1,A2,A3,A4)R,XKOI,XKOH,NF 
COMMON  /UNLDM/  CMI,CMH,SEMAX,EEMAX,EEP,EEI 
COMMON  /GPARM/  RO,CO,GVO,SS,PB 
CHARACTER* 25  NAMl,NAM2,NAM3tNAM4,NAM5,NAMG 
C 

C  READ  INPUT  AND  OUTPUT  FILE  NAMES 
C 

WRITE(*, *)  ’Type  input  file  with  SKELETON  STRAIN-STRESS  pairs’ 
READ(*,’(A) ’)  NAM2 

WRITE(*,*)  ’Type  input  file  with  EFFEC.  VER.  STRAIN,  vs  Ko’ 
READ(*,’(A)’)  NAM4 

WRITE (*,*)  ’Type  output  file  for  (STRN(R) ,SE,PF,S,LF) ’ 

READ(*, ’ (A) ’)  NAM3 

WRITE (*,*)  ’Type  output  file  for  (STRN(R)  ,SEH,SH,BMF,POR)  ’ 
READ(*,’(A)’)  NAM5 

WRITE (*,*)  ’Type  output  file  for  (EV(S) ,P,BMX,PQR) ’ 

READ(*,’(A)’)  NAM8 

OPEN (8 , FILE=NAM3 , STATUS= ’ NEW ’ ) 

OPEN  (3 ,  FHE=NAM5 ,  STATUS*  ’ NEW ’ ) 

OPEN (66 , FILE=NAM6 , STATUS* ’ NEW ’ ) 

C 

C  READ  VARIABLE  Ko  AS  A  FUNCTION  OF  EFFECTIVE  VERTICAL  STRAIN 
C 

OPEN (8 , FILE=NAM4 , STATUS* ’ OLD  * ) 

RBAD(8,*)  NKO 
DO  50  1=1, NKO 
50  READ (8,*)  EKO(I) ,RKO(I) 

CLOSE (8) 

C 

0  READ  NF 
0 
C 

C  NW 

C 

C  EMAX 

0  NDIV 

0 

C  NDIVU 

C  NDIVM 

0  LPRINT 

C 

C  PF 

C  SE 

C  V 

C  POR 

C 
C 

OPEN (55 , FILE=NAM2 , STATUS* ' OLD ’ ) 

C 

READ (55,*)  NF,NW, EMAX, NDIV, NDIVU, NDIVU .LPRINT 

0 

BEAD(55,*)  PF,SE,V,POR 
C 


0  -  PRESSURE  IN  PSI 

1  PRESSURE  IN  BARS 

2  -  PRESSURE  IN  PSI  AND  Kw  =  0.3E6,  Kg  *  5.E6 

0  ---  FRESH  WATER 

1 - SEA  WATER 

MAXIMUM  STRAIN 

NUMBER  OF  DIVISION  BETWEEN  SKELETON  DATA  POINTS 
DURING  LOADING 

NUMBER  OF  DIVISION  FOR  UNLOAD  PRIOR  TO  LIQUEFACTION 
NUMBER  OF  DIVISION  FOR  UNLOAD  AFTER  LIQUEFACTION 
OUTPUT  PRINT  INTERVAL 

INITIAL  PORE  WtfER  PRESSURE 
INITIAL  EFFECTIVE  VERTICAL  STRESS 
INITIAL  POISSON’S  RATIO 
VIRGIN  POROSITY 
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O  O  O  Cl  CIO  OOOOQ  O  O  Cl  Cl  O  O  O  Cl  Cl  Cl  Cl  O  O  Cl  O  O  O  O  O  Cl  Cl  Cl  Cl  O  O  Cl  Cl  Cl 


READ  UNLOADING  CONSTANTS  FOR  SKELETON 

A1,A2  :  PARAMETERS  FOR  RECOVERABLE  STRAINS 
A3,A4  :  PARAMETERS  FOR  INITIAL  UNLOADING  MODULUS 
R  :  FRACTION  OF  STRESSES  W.R.T.  HOOK 
XKOI  :  KO  FOR  INITIAL  PORTION 
XKOH  :  KO  FOR  HOOK  PORTION 


READ (55 , *)  A1,A2,A3,A4,R,XKOI,XKOH 


READ  SOLID 
RO 
CO 
GVO 
SS 
PB 


GRAIN  PROPERTIES 

:  INITIAL  GRAIN  MASS  DENSITY 

:  INITIAL  WAVESPEED 

:  POISSON’S  RATIO  AT  LOW  STRESS 

:  SLOPE  OF  PROPAGATION  VELOCITY  TO  PARTICLE  VELOCITY 

:  THRESHOLD  PRESSURE 


note  that  units  are  lb,  second  and  inch  for  grain  properties 


READ (55,*)  R0,C0,GV0,SS,PB 


READ  EFF.  STRAIN-STRESS  DATA  FOR  SKELETON 


READ (55,*)  NPOINT 
DO  120  1=1, NPOINT 
120  READ (55,*)  RR(I),CM(I) 
CLOSE (55) 

SET  INITIAL  CONDITIONS 


NFIRST  =  0 
LF  =  0 


LF  =  0  :  VIRGIN  LOADING 

1  :  UNLOADING  (INITIAL  PORTION) 

2  :  UNLOADING  (HOOK  PORTION) 

3  :  UNLOADING  (LIQUEFIED) 

XF(SB.EQ.O.O.AND.PF.EQ.O.O)  GO  TO  142 
FOR  NON-ZERO  IN-SITU  STRESSES 
COMPUTE  INITIAL  EFFECTIVE  STRAIN  (EE) 

I  -  1 

130  CONTINUE 

XF(SE.GE.CM(I) .AND.SE.LE.CM(I*1))  GO  TO  140 
I  =  1*1 
GO  TO  130 
140  CONTINUE 

EE  =  RR(I) ♦ (HR (1*1) -RR(I) ) *  (SE-CM(I) ) / (CM(I+1) -CM(I) ) 

KK  »  I 
RR(I)  =  EE 
CU(I)  =  SE 

COMPUTE  INITIAL  WATER  STRAIN  (EWI) ,  ASSUME  ELASTIC 

CALL  BULKW(PF,HKW,NFINW) 

EWI  =  PF/BKW 

COMPUTE  GRAIN  STRAINS  EGI  AND  EG  DUE  TO  GRAIN  PRESSURE (PCI) 
AND  PORE-WATER  PRESSURE (PF) ,  ASSUME  ELASTIC 
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A  O  O  Q  <5000  oo  o  o  o  o  o  o  o  o  o  o  o  o 


PGI  ==  PF+(1.+V)*SE/(3.*(1.-P0R)*(1.-V)) 
CALL  BULKG (PGI , 0 . 0 , POR , BKGI , NF) 

COMPUTE  SECANT  SKELETON  BULK  MODULUS,  BSKN 


BSKN  =  (SE/EE)*(1./3.)*(1.+V)/(1.-V) 

PT  =  (2./3.)*((l.-2. *V) /  (1 . -V) ) * (BSKN /BKGI) *PF 
PGI  =  PGI-PT/(1.-P0R) 

EGI  =  PGI/BKGI 

CALL  BULKG (PF, 0.0, POR, BKGF,NF) 

EG  =  PF/BKGF 

COMPUTE  RELATIVE  VOLUMES 
ASSUME  INITIAL  SAMPLE  VOLUME  =1.0 

VGO  =  1 . -POR 

VGI  =  VGO* (1. -EGI) 

VI  =  1 .  -EE-  (1 .  +V)  *EG/  (3 .  *  (1 .  -V) ) 

VWI  =  Vl-Vfll 

VWO  =  VWI/(1.-EWI) 

VO  =  VGQ+VWG 

COMPUTE  INITIAL  POROSITIES,  STRAINS  AND  STRESSES 

PORI  =  VWO/VO 
POR  =  VWI/VI 

STRN  =  l.-YI 

S  =  PP>SE 

PE  =  (l.*V) *SE/(3.*(1.-V))-PT 
GO  TO  146 

FOR  ZERO-INITIAL  STRESSES 

142  CONTINUE 
VO  =  1.0 
PORI  «  POR 
STRN  =  0.0 
EB  =  0.0 
S  =>  0.0 
SE  =  0.0 
PE  =  0.0 
PF  =  0.0 
KK  =  1 

146  CONTINUE 

STEP-BY-STEP  CALCUUTIONS 


L  a  0 
C 

150  L  =  L+l 

IF(NFIRST.EQ.O)  GO  TO  200 
K  =  K*1 

IF(K.EQ.NDIV)  GO  TO  200 
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o  o  o  ooo  o  o  o  ooo  o  o  o  o  o  o  o  o o  o  o  o  ooo 


GO  TO  300 
200  CONTINUE 
L  =  LPR.INT 
K  =  0 
NFISST  =  1 
KK  =  KK+1 

DDTRN  =  (RR(KK) -RR(KK-1) ) /NDIV 
300  CONTINUE 

CALCULATE  SKELETON  CONSTRAINED  MODULUS 

CALL  CONMS(EE,BMS,LF) 

CALCULATE  SKELETON  BULK  MODULUS 

CALL  BULKS (BMS , BKS , EE , LF , AK) 

CALCULATE  GRAIN  BULK  MODULUS 
PORG=POR 

IP(NF.EQ.2)  PORG=PORI 

CALL  BULKG (PF ,PE , PORG , BKG , NF) 

CALCULATE  WATER  BULK  MODULUS 

CALL  BULKW(PF,BKW,NF,NW) 

CALCULATE  UNDRAINED  BULK  MODULUS 

AA  =  (1 . -PORI) / (l . -POR) 

XF(NF.EQ.2)  AA  =  1.0 

BMF  a  BMS+ (BKG-AA*BKS) / (1 . *PQRI*BKG* (BKG-BKW) / 
.  (BKW* (BKG-AA*BKS) ) ) 

CALCULATE  TOTAL  STRAIN  INCREMENT 

DSTRN  =  DDTRN*(BKG~BKS)/(BKG-BMF*AK) 

CALCULATE  PRESSURE  INCREMENTS 

DS  «  BMP* DSTRN 

DPF  a  DSTRN* (BMP-BMS) /  (1 , -BKS/BKG) 

DSE  =  DS-DPF 

DPE  a  BKS* (DSTRN-DPF/BKG) 

UPDATE  PRESSURES 

PF  a  PF+DPF 
PB  =  PE* DPE 
S  =  S*DS 
SE  a  SE*DSE 

SEH  a  1.6  *  PE-  0.5  •  SE 
SB  =  SEH+PF 

UPDATE  STRAINS  AND  POROSITY 


BE  a  EE*DDTRN 
STRNO  =  STRN 
PORO  =  POR 
STaN  =  STRN*DSTRN 
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DVW  =  PORO* (1 . -STRNO) -PGRI*DPF/BKW 
POR  =  DVW/(1.-STRN) 

C 

IFfLF.EQ.O.AND.STRN.GE.EMAX)  GO  TO  360 
IP(L.NE.LPRINT)  GO  TO  380 
360  L  =  0 

WRITE  TOTAL  VERTICAL  STRESS:  S 

TOTAL  HORIZONTAL  STRESS  (SH) 

UNDRAINED  CONSTRAINED  MODULUS:  BMP 
EFFECTIVE  VERTICAL  STRESS:  SE 
EFFECTIVE  HORIZONTAL  STRESS  (SEH) 

PORE  WATER  PRESSURE:  PF 
POROSITY:  POR 

AS  A  FUNCTION  OF  TOTAL  STRAIN:  STRN 

370  CONTINUE 

TEMP  =  STRN* 100. 

WRITE (6, 2001)  TEMP,SE,PF,S,LF 

2001  FORMAT (4E12 . 4 , 15) 

WRITE (3, 2002)  TEMP .SEH.SH.BMF, POR 

2002  FORMAT (5E12. 4) 

C 

380  CONTINUE 
C 

IFfLP.EQ.O.AND.STRN.GE.EMAX)  GO  TO  400 
IF (LF. GT.O. AND. LF.LT. 3. AND.SEH.LT. 0.0)  GO  TO  440 
IF(LF. EQ.l. AND. K.EQ. (NDIV-1))  GO  TO  420 
XF(LF.EQ.2.AND.K.EQ. (NDIV-1))  GO  TO  440 
IF(LF.EQ.3)  GO  TO  500 
GO  TO  150 
C 

400  CONTINUE 
NDIV  s  NDIVU 
SEMAX  **  SE 
EEMAX  «  EE 
CALL  SKUNL 

0 

IF(R.GT.O.O.AND.R.LT.l.O)  GO  TO  414 
C 

LF  b  2 

RR(KK)  e  EEMAX 
RR(KK+1)  =  EBP 
K  =  NDIV-1 
GO  TO  ISO 
C 

414  CONTINUE 
LF  =  1 

RR(KK)  *  EEMAX 
RR(KK+1)  =>  EBI 
RR(KK*2)  =  EEP 
K  =  NDIV-1 
GO  TO  ISO 
C 

420  CONTINUE 
LF  *  2 
K  =  NDIV-1 
GO  TO  150 
C 
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440  CONTINUE 
LF  =  3 

EVMAX  =  STRN 

CALL  MKTU  (EVMAX ,  NDIYM ,  PORI  ,POR,  BMF ,  PP ,  NF ,  LPRINT ,  NW) 

S  =  PF 

SH  =  PF 

SE  =  0.0 

SEH  =  0.0 

GOTO  370 


C 

500  CONTINUE 
C 

CLOSE  (6) 

C 

CLOSE (6) 

CLOSE (3) 

CLOSE (66) 

STOP 

END 

Q  *******«*»*****************+*+********************«****** 

SUBROUTINE  BULKG (PF , PE , PQR , BKG , NF) 

COMMON  /GPARM/  RO.CO.GVO^.PB 
C 

C  ,  GRAIN  BULK  MODULUS 
C 

IP(NF.NE.2)  GO  TO  50 
BKG  =  5000000. 

RETURN 

G 

50  CONTINUE 

P  «  PF+PE/ (1 . -POR) 

0 

PG  =  P 

IP(NP.Eq.l)  PG  a  P*14,S 

0 

0  COMPUTE  TANGENT  CONSTRAINED  MODULUS,  CG 
C 

SROC  «  SQRT((RO*CQ) **2*4. •BQ*S8*PQ) 

SFP  =  SROC-RO«CO 

?P  =  RO*CQ*CO+CO*SFP+ (SFP*SFP) / (4 . *R0) 

0 

SFPP  =  2.*R0*SS/SR0C 

PPP  a  C0*SFPP*SFP*SFPP/(2.*RQ) 

0 

CG  a  (FP*FP) / (FP-PG*FPP) 

C 

C  COMPUTE  TANGENT  BULK  MODULUS,  BKG 
0 

IP(PG.GT.PB)  GO  TO  100 
C 

GP»(2 . /3 . ) • ( (1 . -2 . *GVO) / (1 . -GVO) ) * (PG/PB) *  (1 . ♦GVO) /  (3 ,  • (1 . -GVO) ) 
C 

GG  TO  200 
C 

100  CONTINUE 
C 

GP  =  1.0 
C 

200  CONTINUE 


C 
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BKG  =  GP*CG 
C 

IF(NF.EQ.l)  BKG  =  BKG/14.5 
C 

RETURN 

END 

C  ********************************************************* 
SUBROUTINE  BULKW(PF,BKW,NF,NW) 

C 

IF(NF.NE.2)  GO  TO  50 
BKW  =  300000. 

RETURN 

C 

50  CONTINUE 
C 

P  =  PF 

IF(NF.EQ.O)  P  =  PF/14.5 
C 

PW1  =  10000. 

PW2  =  100000. 

C 

IF(NW.EQ.l)  GOTO  40 
C 

C  FRESH  WATER  (NW=0) 

C 

Cl  *  0.00031775 
C2  =>  6.2713 
C3  =  21746.9 
C 

C4  =  9.49 
C5  a  21920. 

C 

ce  =  o.oooooeas 

07  =»  1.337 
C8  a  768348. 

GOTO  60 
C 

C  SEA  WATER 

C 

40  CONTINUE 
C 

Cl  »  0.000342831 
C2  «  6.18393 
C3  a  23806.8 
C 

C4  b  9.6546 
CS  a  24758. 

C 

C6  a  0.00000291437 
07  a  3.36508 
08  a  614918. 

C 

60  CONTINUE 

0 

IF(P.GT.Pfl)  GO  TO  100 
C 

BKW  »  C1.P*P*C2*P*C3 
GO  TO  300 

100  IF(P.GT.PW2)  GO  TO  200 
G 
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o  o  o 


BKW  =  C4+P+C5 
GO  TO  300 
C 

POO  CONTINUE 

BKW  a  CS*P*P+C7*P+C8 
C 

300  CONTINUE 
C 

ZFCMF.SQ.l)  RETURN 
C 

BKf  a  BKW*14.5 
C 


C 


C 


RETURN 

END 


*******+****************+**i,mt:mmmttm*********m* 

SUBROUTINE  CQNMS (EE.BMS.LF) 

COMMON  /SKTON/  RR (100), CM (100) 

COifiiON  /UNLDM/  CMI,CMH,SEMAX,EEMAXJEEP.EEI 
IF(LP.NE.O)  GO  TO  300 
1  =  1 

100  CONTINUE 

IF(EE.GE.RR(I).AND.EE.LE,RR(I+l))  GO  TO  200 
I  =  !♦! 


GO  TO  100 
200  CONTINUE 

^(CUd.l)  -CM  (I) )  /  (RR(I+1)  -RR(I) ) 

300  EP(LF.NE.l)  GO  TO  400 
BMS  =  CMI 
RETURN 

400  IP(LP.NB.2)  GO  TO  600 
BUS  a  CUE 
RETURN 

500  CONTINUE 
BUS  »  0.0 
RETURN 
END 


SUBROUTINE  BULKS(BUS,BKS,EE,LP,A) 


FIND  Ko 


IP(LP.EQ,3)  GO  TO  100 
CAJLL  CALKO(EB,XKO,LF) 

A  «  (2.0*)CKOd.Q)/3. 

BKS  a  A* BMS 
RETURN 

400  CONTINUE 
BKS  =  0.0 
RETURN 
END 

,XKO,LF) ****************** 

0) ,RX0(1OQ) 
,A3,A4,R,XK0I,XKQB,NF 
00 

100  CONTINUE 

^(EE.GE.EKOdJ.AND.EE.LE.EKOCI^l))  GO  TO  200 


SUBROUTINE  CALKQ(EB 
COMMON/KODAT/EKO (10 
COMMON/UNLDP/  A1.A2 
IF(LP.NE.O)  GO  TO  3< 
I  =  1 
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GO  TO  100 
200  CONTINUE 
XKO'  =  RKO(I) 

RETURN 

300  IF(LF.NE.l)  GO  TO  400 
XKO  =  XKOI 
RETURN 

400  IP(LF.NE.2)  GO  TO  300 
XKO  =  XKOH 
500  RETURN 
END 

C  ********************************************************* 

SUBROUTINE  SKUNL 

COMMON/UNLDP/  A1 , A2 , A3 , A4 , R , XKOI , XKOH , NF 
COMMON/UNLDU/  CMI,CMfl,SEMAX,EEMAX,EEP,EEI 
C 

C  CALCULATE  UNLOADING  MODULUS 
C 

CMI  =  A3*(SEMAX**A4) 

C 

IF(R.GT.O.O.AND.R.LT.l.O)  GO  TO  100 

0 

EEE  -  SEMAX/CUI 
EEP  =  EEMAX-EEB 
C 

CMH  =  cm 

XKOH  «  XKOI 
EEI  =  EEP 
C 

RETURN 

C 

100  CONTINUE 
C 

C  CALCULATE  RECOVERABLE  (BEE)  AND  PERMANENT  (EEP)  STRAINS 
C 

EEE  a  Al*(SEMAX*»A2) 

EBP  a  EEMAX-EEE 
C 

CMA  =  SEMAX/EEB 

£81  «  EEMAX-fl . -R) *SEMAX/CMI 

CMH  «  R»(RiA/ (l . - (1 . -R)  ♦CMA/ CMI) 

C 

RETURN 

END 

C 

C - * — - 

c 

SUBROUTINE  MmU(EVMA2,NDlVM,POai,POS,BlCC,PPfNP,LPRLNT,N») 
C 

L  -  LPRINT 
EV  =  0.0 
P  =  0.0 

DEV  =  BVMAX/NDIVM 
FOR  =  PORI 
DO  100  I«1 , NDIVM 

CALL  BULKGfP, 0.0, PORI, BKG.NP) 

CALL  BULKW(P,BKW,NP,NW) 

BMX=BKG*BKW/ (BKW*PORI. (BKG-BKW) ) 

DP*BMX*DBV 

EVO=EV 
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PORO=PGR 

EY=EV+DEV 

DVW=PORO* (1-EVO) -PQRI*DP/BKW 
PQR=DYW/ (1-EY) 

P=P+DP 

IFfl.EQ.NDIVM)  GO  TO  40 
IF(L.NE.LPRINT)  GO  TO  50 
40  WRITE (88, 2000)  EV*100,P,BMX,POR 

L  =  0 

SO  L  =  Lt-1 

100  CONTINUE 
2000  FORMAT (4E12. 4) 

PF=P 

RETURN 

END 
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APPENDIX  B 


EXAMPLE  PROBLEM  USING  PROGRAM  NKOCP 
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Table  B.l.  Listings  of  input  files  for  example  problem 


Input  file:  TAPE  55 


2,  0,  0.03,  100,  100,  100,  10 

0.0,  0.0,  0.2,  0.1 

0.0,  0.0,  4000000.,  0.0,  0.0,  0.25,  0.25 
0.0,  0.0,  0.0,  0.0,  0.0 
3 

0.0  ,  0.0 
0.02  ,  20000. 

0.04  ,  40000. 


Input  file:  TAPE  8 


2 

0.0  ,  0.25 

0,04  ,  0.25 
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Table  B.2  Listings  of  output  files  for  example  problem. 


Outpute  file:  TAPE  6 


Total 

axial 

strain  {%) 

0 . 2446E-01 
0.2691E-00 
0.5137E+00 
0.7583E+00 
0.1003E+01 
0 . 1248E+01 
0 . 1492E+01 
0.1737E+01 
0.1981E+01 
0.2226E+01 
0.2471E+01 
0.2715E+01 
G.296QE+01 
0 . 3009E+01 
0 , 3002E+01 
0.293QE+01 
Q.2859E+01 
0.2788E+01 
0.2719E+0'i 
0.2645E+Q1 
0.2573E+01 
0 . 2S52E+01 


Effective 

axial 

stress  (psi ) 

0.2000E+03 
0.220CE+04 
0 . 4200E+04 
0. 8200E+04 
0.8200E+04 
0.1C20E+05 
0.1220E+05 
0. 1420E+05 
G.1820E+05 
0.1820E+05 
0.2020E+05 
0.2220E+0S 
0 . 2420E+05 
0.2480E+0S 
0 . 2435E+05 
0.2189E+05 
0.1943E+05 
0.1697E+05 
0. 1451E+05 
0.1205E+OS 
0.9594E+04 
O.OOOOE+OO 


Pore 

pressure 

(psi) 

0 . 4463E+03 
0 . 4909E+04 
0.9372E+04 
0. 1383E+05 
0.1830E+05 
0.2276E+05 
0 . 2722E+0S 
0.3169E+05 
0.381SE+0S 
0.4061E+05 
0.4507E+0S 
0.49S4E+05 
0 . 5400E+05 
0 . 5489E+05 
0 . 5479E+0S 
0.5381E+05 
0 . 5282B+05 
0.5183E+05 
0 . S084E+06 
0.4985E+05 
0.4886E+0S 
0.4972E+0S 


Total 

axial 

stress  (psi) 

0.6483E+03 
0.7109E+04 
0.1357E+05 
0 . 2003E+05 
0.2650E+05 
0.3296E+05 
0.3942E+05 
0.4589E+05 
0.5235E+05 
Q.5881E+05 
0 . 6527E+05 
0.7174E+05 
0.7820E+05 
0 . 7949E+0S 
0.7915E+05 
0.7S70E+05 
0 . 7225E-*Q5 
0 . 6880E+0S 
0.6S3SB+05 
0.6191E+CS 
0.5846E+05 
0.4972E+05 


Loading 

Status 

o 

o 

o 

o 

o 

0 

0 

0 

0 

0 

0 

0 

Q 

0 

2 

2 

2 

2 

2 

2 

2 

3 
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Table  B.2.  Listings  of  output  files  for  example  problem  (concluded). 


Output  file:  TAPE  3 

Total  Effective  Total  Undrained 

axial  radial  radial  Modulus 

strain  (%)  stress  (psi)  stress  (psi )  (psi) 


0.2448E-01 
0.2691E+00 
0.5137E+00 
0.7S83E+00 
0.1003E+01 
0.1248E+01 
0.1492E+01 
0.1737E+01 
0.1981E+01 
0.2226E+01 
0.2471E+01 
0.2715E+0X 
0.2960E+01 
0 . 3009E+01 
0.3002E+01 
0.2930E+01 
0.2859E+01 
0 . 2788E+01 
0.2/ J.6E+01 
0.2645E+01 
0.2573E+01 
0.2SS2E+0j. 


0.1653E+02 
0.1818E+03 
0.3471E+03 
0.5124E+03 
0.8777E+03 
0.8430E+03 
0.1008E+04 
0. 1174E+04 
0.1339E+04 
0.1504E+04 
0.1669E+04 
0.183SE+04 
0.2000E+04 
0 . 2033E+04 
0 . 2001E+04 
0.1683E+04 
0 . 1384E+04 
0 . 1046E+04 
0.7273E+03 
0.4088E+03 
0 . 9032E+02 
O.QQOOE+QQ 


0.4628E+03 
0.5091E+04 
0.9719E+04 
0.143SE+05 
0.1898E+05 
0.2360E+05 
0.2823E+05 
0.3286E+0S 
0.3749E+05 
0.4212E+05 
0.4874E+05 
0.S137E+05 
0.5800E+05 
O.S693E+OS 
0.5679E+05 
n . SS49E+0S 
U.5418E+05 
0.5287E+0S 
0.S1S7E+0S 
0 . 5026E+05 
0.4895E+05 
0 . 4972E+Q5 


0 . 2642E+07 
0.2642E+07 
0.2642E+07 
0.2842E+07 
0.2642E+07 
0.2642E+07 
0.2642E+07 
0.2642E+07 
0.2642E+07 
0.2642E+07 
0.2842E+0? 
0.2842E+07 
0 . 2842E+07 
0.2642E+07 
0 . 4831E+07 
0.4831E+07 
0.4831E+07 
0.4831E+07 
0.4831E+07 
0.4831E+07 
0.4831E+07 
0.1948E+07 


Current 

Porosity 


0.9988E-01 
0.9883E-01 
0.9738E-01 
0.9612E-01 
0.948SE-01 
0.9358E-01 
0.9230E-G1 
0.9102E-01 
Q.8S73E-01 
0 . 8843E-01 
0.8713E-01 
0 . 8582E-01 
0.8450E-01 
0.8424E-01 
0.8427E-01 
0 . 8454E-01 
0.84-32E-01 
0.8510E-01 
0.8537E-01 
Q.8565E-0I 
0.8S92E-01 
0.8S61E-01 


Output  file:  TAPE  66  (Mixture  Response) 

Total  Pressure  Undrained  Current 

strain  (%)  (psi)  Modulus  Porosity 


0.2552E-01 

0.2807E+00 

0.5359E+G0 

0.7911F.+00 

0.1048E+01 

Q.13G2E+01 

0.1557E+01 

0.1812E+01 

0.2067E+01 

0.2322E+01 

0.2S52E+01 


0.4972E+03 
0.54696+04 
0.1044E+05 
0. 1541E+05 
Q.2G38E+05 
0.2S35B+05 
0.3033E+05 
0.3530E+05 
0 . 4027E+05 
0.4S24E+0S 
0.4972E+05 


0.1948E+07 
0.1948E+07 
0  ’9488+07 
0. 1U48E+07 
0.1948E+07 
0.  x948E+07 
0.1948E+07 
0.1948B+07 
0.1948E+07 
0. 1948E+07 
0.1948E+07 


0.8985E-01 
0 , 9845E-01 
0.9704E -01 
0.3562E-01 
0.9419E-01 
0.9276B-01 
0.9131E-01 
0.8986E-01 
0.8340E-01 
0.8694E-01 
0.8561E-01 
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Stress 


UNDRAINED  UNIAXIAL  STRAIN  LOAD-UNLOAD  RESPONSE 
OF  IDEALIZED  SATURATED  POROUS  MATERIAL 


TOTAL  AXIAL  STRESS 


MIXTURE  STRESS 


/  ^ 

y  //' 

/  / /r- PORE  WATER 

/  ^  PRESSURE 

/  // 

/  // 

// 

// 

'/  EFFECTIVE  AXIAL  STRESS 


EFFECTIVE  RAOIAL  STRESS 


Strain  (%) 


Figure  B.l.  Example  problem  using  program  NKOCP. 
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APPENDIX  C 

LISTING  OF  PROGRAM  TWAVE 


PROGRAM  TWAVE 
C 

IMPLICIT  DOUBLE  PRECISION  (A  -  Z) 

PROGRAM  TWAVE  CALCULATES  WAVE  SPEEDS  AND  DAMPING  CONSTANTS 
OF  WAVES  OF  FIRST  AMD  SECOND  KIND  IN  THE  POROUS  MEDIUM 
SATURATED  BY  FLUID 

FOR  LOW  FREQUENCY-PERMEABILITY  PRODUCT  (  KAPPA  <  2.0  ) 

USE  :  RB  o  0.0 

FOR  HIGH  FREQUENCY-PERMEABILITY  PRODUCT  (  KAPPA  >  2.0  ) 

USE  :  R  =  0.0 

NF  =  0  :  EXCITATION  FREQUENCY  IS  GIVEN 

1  :  COMPUTE  APPROXIMATE  CRITICAL  FREQUENCY  AND  USE  IT 

openCunitaS.fi.len'IWAVE'  ,status='old’  ) 
open(unit=6,file=*0WAVEf ,status='new' ) 

CALL  RDATA 

CALL  WAVE 

CALL  OUTPT 

STOP 
END 

SUBROUTINE  RDATA 

IMPLICIT  DOUBLE  PRECISION  (A  -  Z) 

INTEGER  NF 

COMMON  /INP/  W , POR, PER , R , RB , GW ,G , SGF , SGG , BF , BG , BS ,CS , NF 

READ  INPUT  DATA 

READ  (5,1000)  NF 
READ  (5,2000)  W 
READ  (5,2000)  POR,PER,R,RB 
READ  (5,2000)  GW.C.SGF.SGG 
READ  (5,2000)  BF,BG,BS,CS 
C 

1000  F0RMAT(I5) 

C 

2000  FORMAT( 4F10.0 ) 

C 

RETURN 

END 

C 

SUBROUTINE  WAVE 

IMPLICIT  DOUBLE  PRECISION  (A  -  Z) 

INTEGER  NF 

COMMON  /INP/  W , POR , PER , R , RB , GW ,G , SGF ,SGG , BF , BG , BS , CS , NF 
COMMON  / OUT/  CP , CM , ALPHAP , ALPHAM , RK 
C 

ROF  -  SGF*GW/G 
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ROG  =  SGG*GW/G 
RO  =  P0R*R0F+( 1 .-P0R)*R0G 
C 

A2  =  R0F*(1.+R)/P0R 

A3  =  P0R/BF+(1.-P0R)/BG-BS/(BG*BG) 

A4  =  l.-BS/BG 
C 

B1  =  R0*A2*A3-(R0F**2)*A3 
B3  =  2.*ROF*A4-A2*(A4**2)-A2*A3*CS-RO 
C 

IF(NF.EQ.O)  GO  TO  50 

COMPUTE  CRITICAL  FREQUENCY 

B40  =  -(A4*A4+A3*CS) 

El  =  (B3/CS)**2-4.*(B1/CS) 

E2  =  (B40/CS)**2 

WC  .  (GW/PER)*DSQRT(E2/E1 ) 

W  =  WC 

50  CONTINUE 

COMPUTE  WAVE  SPEEDS  AND  DAMPING  CONSTANTS 
ONE  =  1.0 

PIE  o  4.0*DATAN(0NE) 

A 1  -  -GW/(PER*W) 

B2  -  R0*A 1*A3 
B4  »  Al*( A4*#2)+A1*A3*CS 

FRICTION  CORRECTION  FOR  HIGH  FREQUENCY  MODE  (  RB  .NE.  0.0  ) 
IF(RB.EQ.O.O)  GO  TO  80 
RKO  «  DSQRT(8.*W*PER/(P0R*G)) 

RK  »  RB*RK0 

IF(RK.LT.20. )  CALL  FKEXC( RK , FI , F2) 

IF(RK.GE.20.)  CALL  FKAPP( RK , FI , F2 ) 

CALCULATE  NEW  B1,B2,B3»AND  B4 

B2B  «  B2 
BAB  «  B4 
C 

B 1  »  B1-B2B*F2 
B2  »  B2B*P1 
C 

B3  a  B3+B4B>!,F2 
B4  a  B4<*FI 
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C 


80  CONTINUE 


C 

C 


c 

85 


C 

90 

C 


C 


Cl  =  (B3**2-B4**2)/(CS**2)-4.*B1/CS 
C2  =  -4.*B2/CS-2.*B3*B4/(CS**2) 


DAPZERO  =  1.0D-100 
C1ABS  =  DABS (Cl ) 

IF(C1ABS.GE. DAPZERO)  GO  TO  85 
IF(C2.GE.0.0)  T2  =  0.5*PIE 
IF(C2.LT.0.0)  T2  =  1.5*PIE 
GO  TO  90 


CONTINUE 

C12  =  DABS(C2/C1 ) 

AT2  =  DATAN(C12) 
IF(C1.GE.0.0.AND.C2.GE.0.0)  T2 
IF(C1.LE.O.O.AND.C2.GE.O.O)  T2 
IF(C1.LE.0.0.AND.C2.LE.0.0)  T2 
IF(C1.GE.0,0.AND.C2.LE.0.0)  T2 


AT2 

PIE-AT2 

PIE+AT2 

-AT2 


CONTINUE 


T12  »  0 . 5*T2 

Z4  «  (C1**2+C2J!'*2)**0.25 

DIP  -  -B3/CS+Z4*DC0S (T1 2 ) 
DIM  -  -B3/CS-Z4*DC0S(T12) 
D2P  -  B4/CS+Z4*DSIN(T12) 
D2M  »  B4/CS-Z4*DSIN(T12) 
D3  »  2 .  *B  1  /CS 
D4  i  2 . *B2/CS 


COMPUTE  WAVE-SPEED 


CALL  WDAMP( DIP , D2P, D3 , D4 , W ,CP , ALPHAP) 

C 

CALL  WDAMP( DIM , D2M , D3 , D4 ,W ,CM , ALPHAM) 

C 

Cl  -  CP 
A 1  »  ALPHAP 

IF(CP.GE.CM)  GO  TO  100 
CP  «  CM 

ALPHAP  «  ALPHAM 
CM  =>  Ci 
ALPHAM  a  A1 
100  CONTINUE 
C 

RETURN 

END 

C 

SUBROUTINE  OUTPT 

IMPLICIT  DOUBLE  PRECISION  (A  -  Z) 

INTEGER  NF 

COMMON  /I NP/  W , POR, PER , R, RB,GW,G , SGF.SGG , BF , BG , BS ,CS,N F 
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cc:mon  /OUT/  CP, CM, ALPHAP,ALPHAM,RK 
PRINT  INPUT  DATA  WITH  LABELS 

WRITE( 6 , 2000)  NF , W, POR , PER , R , RB , GW, G , SGF , SGG , BF, BG , BS , CS 
FORMAT  STATEMENT 


2000  F0RMAT(  // 

.  '  OPTIONAL  FLAG - -  *  ,15  / 

.  ’  . EQ . 0  :  COMPUTE  AT  GIVEN  FREQUENCY  '  / 

.  '  . EQ . 1  :  COMPUTE  AT  CRITICAL  FREQUENCY  '  // 

.  *  FREQUENCY - -  *  ,  El  3 . 6  // 

.  '  POROSITY - -  ’  ,  FI  2 . 5  // 

.  *  COEFFICIENT  OF  PERMEABILITY - =  \E13.6  // 

.  ’  MASS  INCREMENT  FACTOR - -  \F12.7  // 

.  '  PORE  SHAPE  FACTOR - -  '  ,F12.7  // 

.  '  UNIT  WEIGHT  OF  STANDARD  WATER - -  \E13.6  // 

.  '  GRAVITATIONAL  ACCELERATION - -  \E13.6  // 

.  '  SPECIFIC  GRAVITY  OF  PORE  FLUID - -  *  ,Fi2.5  // 

.  ’  SPECIFIC  GRAVITY  OF  SOLID  GRAIN - -  * , FI 2 . 5  // 

.  '  BULK  MODULUS  OF  PORE  FLUID - -  \E13.6  // 

.  '  BULK  MODULUS  OF  SOLID  GRAIN - -  \E13.6  // 

.  '  BULK  MODULUS  OF  SOLID  SKELETON - =  *  ,  El 3 . 6  // 

.  '  CONSTRAINED  MODULUS  OF  SOLID  SKELETON - =  * , El 3 . 6  //) 


WRITE  COMPUTED  WAVE-SPEEDS 


WRITE( 6 , 2100)  CP.CM.ALPHAP, ALPHAM.RK 


2100  FORMAT ( 

// 

.  '  SPEED  OF  FIRST  KIND  WAVE  — - - 

a 

' ,E13.6 

// 

.  '  SPEED  OF  SECOND  KIND  WAVE  - - - 

—  —  —  3 

' ,E13.6 

// 

.  '  DMAPING  CONSTANT  OF  FIRST  KIND  WAVE  - 

—  —  —  a 

' ,E13.6 

// 

.  '  DAMPING  CONSTANT  OF  SECOND  KIND  WAVE  — 

—  —  —  C£3 

' . E13.6 

// 

.  ’  NON-DIMENSIONAL  PARAMETER,  KAPPA  . 

—  —  —  e 

\EI3.6 

//) 

RETURN 

END 


SUBROUTINE  WD AMP (D1,D2,D3,D4,W,C, ALPHA) 
IMPLICIT  DOUBLE  PRECISION  (A  -  Z) 

C 

D34  .  D3^D3+D4^D4 
C 

K3  =  (Dl*D3-t-D2*D4)/D34 
C 

K4  =  (D2<‘D3-D1‘!<D4)/D34 
C 


C 


SQK34  =  DSQRT(K3*K3+K4!!‘K4) 
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C 

K 1  =  (l./SQRT(2.))*DSQRT(SQK34+K3) 

C 

K2  =  -(l./SQRT(2.))*DSQRT(SQK34-K3) 


C  =  K1+K2*K2/K1 
C 

ALPHA  =  K2*W/(C*K1 ) 

C 

C 

RETURN 

END 

SUBROUTINE  FKEXC(RK , FI , F2) 
IMPLICIT  DOUBLE  PRECISION  (A  -  Z) 
C 

CALL  BESSL(RK,X1,X2,X3,X4) 

C 

CALL  C0MPX(X1,X2,X3,X4,T1,T2) 

C 

FIB  =  1.-2. *T2/RK 
F2B  =  2.*T1/RK 
F33  =  0. 25*RK*T1 
F4B  =  0 . 25*RK*T2 
C 

CALL  C0MPX( F1B,F2B,F3B,F4B,F1,F2) 
C 

RETURN 

END 

C 

SUBROUTINE  BESSL( RK ,X1,X2,X3,X4) 
IMPLICIT  DOUBLE  PRECISION  (A  -  Z) 
INTEGER  NE , N 
C 

NE  «  60 
C 

XI  »  0.0 
X2  a  0.0 
X3  a  0.0 
X4  a  0.0 
C 

FACEV  .  1.0 
N  a  0 
SN  =  -1.0 
SIGN  a  l.o 
C 

100  CONTINUE 
C 

N  a  N+l 

c 

FACEV  .  ( 2 . *N-2 .  J^FACEV 
IF(N.EQ.l)  FACEV  a  1.0 
FACOD  a  (2.*N-1.)*FACEV 
C 

PI  »  (0.5*RK)**(4.*N-4.) 
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P2  =  (0.5*RK)**(4.*N-2.) 

P3  =  (0.5*RK)**(4.*N-5.) 

P4  =  (0.5*RK)**(4.*N-3. ) 

C 

XI  =  X1+SIGN*(P1/FACEV)/FACEV 
X2  =  X2+SIGN*( P2/FAC0D) /FACOD 
C 

X  3  =  X3+SIGN*(P3/FACEV)*(2.*N-2.)/FACEV 
X4  =  X4+SIGN*(P4/FAC0D)*(2.*N-1.)/FAC0D 
C 

SIGN  =  SIGN*SN 
FACEV  =  FACOD 
C 

IF(N.LT.NE)  GO  TO  100 
C 

RETURN 

END 

C 

SUBROUTINE  COMPX ( El , E2 , E3 , E4 , Y 1 , Y 2 ) 
IMPLICIT  DOUBLE  PRECISION  (A  -  Z) 

C 

DD  «  E1*E1+E2*E2 
C 

Y 1  a  (E3*E1+E2*E4)/DD 
C 

Y2  a  (E1*E4-E2*E3)/DD 
C 

RETURN 

END 

C 

SUBROUTINE  FKAPP(RK,F1 ,F2) 

IMPLICIT  DOUBLE  PRECISION  (A  -  Z) 

C 

FA  »  2.9213 

FB  =  1 ,/(4.*SQRT(2.)) 

C 

FBK  -  (FB*RK)**2 
FI  =  DSQRT( FA+FBK ) 

C 

F2  •  FB*RK 
C 

RETURN 

END 


APPENDIX  0 


CALCULATION  OF  EQUIVALENT  COEFFICIENT  OF  PERMEABILITY  FOR  A  FLAT  DUCT 

Figure  8.10b  shows  subdivisions  of  the  cross  section  used  to  compute  the 
equivalent  coefficient  of  permeability  for  the  flat  duct.  The  total  cross 
section  is  assumed  to  be  divided  into  the  middle  section  and  two  square  side 
sections  connecting  the  parallel  faces  of  the  middle  section. 

For  the  middle  section,  we  assumed  the  flow  condition  as  in  a  theoretical 
fist  duct  having  infinite  width.  For  the  side  sections,  the  flow  conditions 
are  assumed  to  be  governed  by  the  equivalent  circular  duct  having  the  same 
perimeter  as  the  three  wet  sides  of  the  square. 

Assuming  that  the  flow  velocity  in  the  middle  section  (area  Aj)  is 
Vj  and  the  flow  velocity  in  each  side  section  (area  A2)  is  V2,  then  the 
total  flow  is  given  by 

Q  •  AjVj  ♦  2  A2V2  (0-1) 

Velocities  are  related  to  the  hydraulic  gradient,  i,  as 

Vj  •  i  (0-2) 


and 


V2  ■  k2  i 


(0-3) 


where  kj  and  k2  are  the  coefficients  of  permeability  for  the  middle  and  the 
side  sections,  respectively,  kj  is  computed  from  the  theoretical  flat  duct 
equation  as 


kl 


»12  Yf 

V 


(0—4 ) 
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where 


=  half  the  height  of  the  duct. 


k2  is  computed  from  the  circular  duct  flow  (Equation  8.6)  with  the  radius  of 
the  duct,  a,  given  by 


a  » 

From  Figure  8.10a,  the  total  flow  is 

Q 

where  A  =  Aj  ♦  2  Aj  (0-7 ) 

and  V  *  keqi  (0-Q) 

keq  represents  the  equivalent  coefficient  of  permeability  for  the  total  cross 
section  of  the  flat  duct.  Substituting  Equation  0-2  and  0-3  into  0-t, 

Q  •  (Ajkj  ♦  2  A2k2)i  (0-9) 

and  substitution  of  Equation  D-7  and  0-8  into  0-6  yields 


3  JT2 
2  ~W~ 


coraputed  by 


=  AV 


(0-5) 


(0-6) 


Q  *»  (Aj  ♦  2A2)  keqi  (0-10) 

From  Equations  0-9  and  0-10,  the  equivalent  coefficient  of  permeability  is 
computed  as 


(0-11) 
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